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\  Abstract 

*• 

The  theory  of  high-gain  error  actuated  feedback  control,  developed 
by  Porter  and  Bradshaw,  is  applied  to  the  design  of  various  lateral-direc¬ 
tional  decoupling  flight  control  systems  for  an  advanced  aircraft.  The 
controllers  developed  in  this  report  utilize  output  feedback  with  propor¬ 
tional  plus  integral  control  to  produce  desirable  closed-loop  responses 
with  minimal  interactions  between  outputs.  Because  of  the  structure  of  the 
system,  measurement  variables  in  addition  to  the  outputs  are  necessary  to 
apply  this  method. The  report  examines  controller  design  robustness  by  varying 
the  flight  conditions  or  maneuver  commands  from  the  ones  the  controller  is 
specifically  designed  for,  and  then  judges  system  performance.  The  resiilts 
show  that  the  controller  is  robust  with  respect  to  varying  flight  conditions, 
but  is  not  robust  with  respect  to  varying  maneuver  commands.  This  report 
also  examines  the  effect  of  first-order  actuator  dynamics  in  the  system  model. 
Actuator  dynamics  are  shown  to  significantly  affect  the  control  system  response 
indicating  that  a  simplified  model,  without  actuators,  is  not  desireable  in 
one's,  control  design  scheme.  Also  a  computer  nrogram  to  determine  transmission 
zeros  and  decoupling  zeros  is  developed. 


HIGH-GAIN  ERROR  ACTUATED  FLIGHT 
CONTROL  SYSTEMS  FOR  CONTINUOUS  LINEAR 
MULTIVARIABLE  PLANTS 

I.  Introduction 

Advamced  Fighter  Technology  Integration  (AFTI)  is  the  culmination  of 
the  most  recent  technology  advemces  in  flight  control  emd  avionics.  The 
AFTI-16  aircraft,  which  incorporates  this  technology,  is  presently  xander- 
going  flight  tests.  A  new  digital  flight  control  system,  along  with 
additional  control  surfaces,  such  as  vertical  ceuinaurds,  allow  the  aircraft 
to  perform  direct  force,  decoupled  maneuvers. 

This  report  contains  an  investigation  of  a  new  controller  design 
method  for  the  AFTI-16.  Specifically,  am  analogue  controller  is  designed 
to  perform  four  lateral-directional  decoupled  mameuvers.  'Riese  mameuvers 
aure  horizontal  translation,  wings  level  tvrn,  constamt  altitude  coordinated 
turn,  and  yaw  pointing.  The  mameuvers  can  be  described  by  three  output 
variables:  side-slip  angle,  3;  roll  amgle,  (j>;  amd  yaw  rate,  r.  Initially 
ad.1  output  variadsles  are  zero.  For  horizontal  tramslation  the  side-slip 
amgle,  B,  chamges,  while  and  r  remain  unchamged.  In  a  wings  level  turn 
the  yaw  rate,  r,  chamges,  while  3  snd  remaLin  unchamged.  In  a  constant 
altitude  coordinated  turn  the  yaw  rate,  r,  is  proportional  to  the  roll 
amgle ,  , 

r  =  ^  sin  (j),  (1-1) 

while  3  remains  zero.  For  yaw  pointing  the  time  rate  of  chamge  of  the 
side-slip  amgle,  3r  must  be  equal  in  magnitude  and  opposite  in  sign  to 


I  ,• 

■ 


the  yaw  rate,  r,  vrtiile  (j)  renuiins  zero. 


Background  , 

Controllers  have  been  designed  for  single  input/single  output  (SISO) 
systems  using  conventional  techniques  for  a  long  time.  The  advance  of 
conputational  technology  readily  permits  the  design  of  controllers  for 
multiple  input/multiple  output  (MXHO)  systems.  A  MIMO  system  design 
offers  a  greater  2unount  of  control  emd  greater  precision.  Conventional 
control  en^loys  widely  accepted  design  techniques,  such  ^is  root- locus 
emd  frequency  response  €malysis.  Multi vauriable  control  design  is  still 
growing  and  a  number  of  different  theories  emd  design  methods  ^ave  been 
developed.  In  this  report  the  design  method  developed  and  proposed  by 
Or  Bricui  Porter  of  the  University  of  Salford,  Englemd  is  utilized  euid 
evaluated. 

Other  multivariable  control  methods  presently  available  utilize  optimal 
control,  entire  eigenstructure  assignment,  and  frequency  domain  techniques. 

In  the  optimal  control  design  method  (Ref. 2)  a  control  law  is  found 
by  minimizing  the  value  of  the  performemce  index  or  cost  functional 


J  »  1/2 


T 

+  u  (t) 


^  ^[y(t)  -  v(t)]^  2^^)  ~  v(t)] 

(t)  R(t)  u(t)  J  dt 


+  l/2[y;(tf)  -  v(tf)]  H  [y.(tf)  -  v(tf)] 
for  the  system  described  by  the  state  emd  output  equations 


(1-2) 


x(t)  »  A  x(t)  +  B  u(t) 


X.(t)  -  C  x(t) 


(1-3) 
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where  y(t)  is  the  output  vector. 

v(t)  is  the  command  vector,  and  2^(t) ,  R(t) ,  auid  H  are  symmetric 
positive  definite  quadratic  weighting  matrices.  The  optimal  control  law 
that  minimizes  J  is 

u*(t)  =  K(t)  x(t)  -  R;^(t)  d(t)  (1-4) 

The  time  varying  feedback  gain  matrix  is 

K(t)  =  -  R"^(t)  P(t)  (1-5) 

where  P^(t)  is  the  solution  of  a  nonlinear  time  varying  matrix  differential 
equation,  commonly  known  eis  the  Riccati  Equation: 

P(t)  =  P(t)  -  P(t)A  +  P(t)B  R“^(t)  B]^P(t) 

-Si  2(t)  (1-6) 

where  P (tf )  =  H 

The  vector  d(t)  is  the  solution  to  the  matrix  differential  equation 
d(t)  =  A^(t)  +  K(t)  B  R"^(t)B^d(t) 

+  2(t)  C,  (1-7) 

\diere  d  (tf )  =  -C_'^  H  ^  v  (tf ) . 

Due  to  the  boundary  conditions  on  the  two  matrix  differential  equa¬ 
tions,  they  must  be  integrated  backwards  in  time.  This  requires  that 
commamd  input  information  be  known  in  the  future.  The  computation  re¬ 
quired  for  the  numerical  solutions  to  these  equations  are  extremely  costly 


euld  difficult  to  achieve.  Also,  selection  of  the  positive  definite 
weighting  matrices,  2.(t) ,  R(t) ,  euid  H,  must  be  made.  This  is  usually 
accon^lished  by  a  trial  zuid  error  iterative  process,  requiring  repeated 
integrations  of  the  matrix  differential  equations. 

In  the  entire  eigenstructxire  etssignment  design  method  (Ref  11),  the 
designer  selects  the  system  closed- loop  eigenvalues,!^,  euid  corresponding 
eigenvectors,  Given  the  system  described  by  eq.  1.3,  a  state  feed¬ 

back  control  law 

u(t)  »  K  x(t)  +J(t)  (1-8) 

is  assumed.  Substituting  this  control  law  into  the  state  equations  yields 
the  closed-loop  system 

x(t)  =  [A  +  B  k]  x(t)  +  B  v(t)  (1-9) 

llie  feedback  ged.n  matrix,  K,  is 

-1 

K  “  f  t  oj’2 »  •  •  • »  iiJjj  ]  •  S2 '  *  •  • '  .§n  ^  (I-IO ) 

The  feedback  geun  matrix  of  eq.  1-10  produces  a  closed-loop  system  with 
the  selected  eigenvalues  and  corresponding  eigenvectors,  provided  all 
eigenvectors,  and  vectors,  satisfy  the  relation 

A.  I  ,  bI  (1-11) 

1  n  ’  — J 

Extreme  care  m\ist  be  tsdcen  vrtien  selecting  the  eigenvectors  from  the  null 

-1 

space  of  I ,  b]  ,  in  order  that  the  matrix  [ ^ »  £2 '  * ’ • 


k 


exists.  Also,  the  eigenvalues  and  eigenvectors  must  be  chosen  to  minimize 
the  amount  of  Interaction  between  outputs.  Selection  of  the  eigenvalues, 
and  eigenvectors,  is  em  art,  requiring  insight  euid  a  lot  of  trial 
2uid  error  iterations. 

Porter's  design  method,  described  in  Chapter  II,  often  produces  a 
successful  controller  design  with  less  trial  amd  error  iterations  than  the 
two  methods  described  cdsove. 

Problem  Statement  amd  Approach 

It  is  desired  to  design  eUid  evaluate  a  flight  controller,  using 
Porter's  design  method  (Ref  7),  for  em  AFTI-16  ad.rcraft  which  simulates 
meuieuvers  of  horizontal  translation,  wings  level  turn,  constemt  altitude 
coordinated  turn,  wd  yaw  pointing.  Ilie  evaluation  consists  of  examining 
robustness  and  effect  of  first-order  actuators  in  the  simulation. 

A  robust  controller  design  C2m  perform  satisfactorily  at  conditions 
other  th2ui  the  one  it  is  specif ic€illy  designed  for.  In  this  report 
robustness  is  exeunlned  by  using  a  controller  designed  for  a  yaw  pointing 
maneuver  at  a  specific  design  flight  condition,  and  using  it  without  chemge 
in  the  simulation  of  a  yaw  pointing  meuieuver  at  other  flight  conditions. 
Also  the  same  controller  is  used  in  the  simulation  of  the  three  other 
m^meuvers  at  the  design  flight  condition.  The  performance  of  the  three 
lumeuvers  is  conpeured  with  the  perfomumce  of  controllers  designed  to 
perform  the  specific  nuuieuver. 

Three  different  controllers  are  designed  to  perform  the  yaw  pointing 
maneuver.  One  at  a  low  dynamic  pressure,  Q,  design  flight  condition, 
wother  at  a  medium  dynamic  pressure  design  flight  condition,  emd  the 


third  at  a  hig^  (^namic  pressure  design  flight  condition.  Each  of  these 
controllers  are  then  included  in  the  simulation  of  a  yaw  pointing  mcuieuver 
at  nine  other  flight  conditions,  covering  a  spectrum  of  dynamic  pressures. 

A  successful  design  is  based  on  several  criteria.  In  this  report  the 
following  specifications  are  established:  First,  the  closed-loop  system 
must  be  staOile.  Second,  all  tr^ulsient  responses  should  be  within  98%  of 
their  steady-state  value  in  three  seconds.  "Wiird,  the  steacty-state 
response  should  achieve  the  commanded  value  with  no  observable  error. 
Fourth,  the  control  surfaces  are  not  commeuided  beyond  a  pre-set  limit. 


Assunmtions 


The  following  assunptions  are  made  to  simplify  the  computations  and 


analysis ; 


1.  The  earth  is  flat  and  non- rotating. 

2.  The  atmosphere  is  at  rest  with  respect  to  the  earth. 

3.  The  acceleration  due  to  gravity  is  constant. 

4.  The  aircraft  is  a  rigid  body. 

5.  The  aircraft  has  a  constant  mass. 

6.  The  2d.rcraft  performs  mameuvers  while  flying  straight  with  wings 
level,  unless  otherwise  specified. 

7.  The  motion  of  the  aircraft  edjout  a  straight  and  level  trim  con¬ 
dition  cem  be  adequately  modeled  by  a  linear  dyneunical  system  of  equations. 

8.  The  dyneunics  of  the  actuators  can  be  adequately  modeled  by  a 
first-order  differential  equation. 

9.  The  side  slip  angle,  ^  ,  cem  be  meeisured. 
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Organization 


First,  the  theory  used  in  porter's  design  is  described,  and  the  method 
for  developing  a  controller  is  discussed.  Also  included  in  Chapter  II  are 
methods  for  including  first-order  actuator  dynamics  in  the  model.  Chapter 
III  presents  the  results  of  applying  this  method  to  an  aircraft  model. 

The  effects  of  varying  flight  conditions  and  memeuver  commcuids  are  described 
in  detail.  Also  the  effects  of  actuator  dynaunics  in  the  simulation  aure 
disciissed.  Chapter  IV  discusses  the  results  cuid  provides  recommendations 
for  further  study. 

Included  in  the  Appendices  are  discussion  of  the  2d.rcraft  equations 
of  motion,  a  computer  program  to  calculate  system  zeros,  amd  the  computer 
program  used  to  perform  design  iterations  and  simulations. 
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Theoretical  Development 


Porter's  design  methods  (Ref  7)  are  both  conceptually  and  computa¬ 
tionally  simple.  These  techniques  utilize  singular  perturbation  methods, 
and  are  applicable  to  the  design  of  both  analogue  and  digital  controllers. 
As  stated  in  the  introduction  these  methods  are  applied  to  tracking  systems 
incorporating  high-gain  error  actuated  analogue  controllers. 

High-gain  frequently  has  the  advantage  of  making  the  closed-loop 
System  robust,  or  insensitive  to  plant  parameter  variations.  Also,  high- 
gain  can  produce  a  tight  tracking  loop  which  may  be  desirable  from  the 
viewpoint  of  disturbance  rejection.  Robustness  is  examined  in  this  re¬ 
port,  but  disturbance  rejection  is  not. 

Error  actuated  control  minimizes  the  difference  between  the  commands 
and  outputs.  The  closed-loop  behavior  becomes  increasingly  tight  and 
non-interacting  as  the  gain  parameter  is  increased. 

Singular  perturbation  methods  are  used  to  exhibit  the  asymptotic 
structure  of  the  transfer  function  matrices.  These  methods  examine  the 
overall  effect  of  letting  certain  parameters  become  smaller  and  smaller. 

By  letting  the  perturbation  parameter,  e,  be  the  reciprocal  of  the  for¬ 
ward  path  gain,  g,  the  results  of  singular  perturbation  theory  can  be 
applied,  since  as  e  — ♦  0,  g  =  1/e—*  “.  With  singular  perturbation  theory 
the  asymptotic  properties  of  the  closed-loop  transfer  functions  can  be 
excunined  to  determine  the  location  of  the  closed-loop  poles,  as  the  gain 
becomes  increasingly  large.  These  asymptotic  closed- loop  poles  determine 


the  system  response. 


The  tracking  system  consists  of  a  linear  multivariedale  plant 
governed  on  the  continuous  time  set  T=  (0,oo)  represented  by  st^te  and 
output  equations  of  the  form  * 


X  (t)  *  A  x(t)  +  B  u(t) 


(t)  =  C  x(t) 


(2.1) 


The  state  and  output  equations,  eq.  2.1,  expressed  in  expanded  form  are 


[^(t)|  r  All  Ai2  xi(t)  0 

X2(t)J  ^2  l2. 


(2.2) 


I  <t)  = 


[£i  £2]  r 


3^(t) 

X2(t) 


(2.3) 


Where  is  a  square  non-singular  matrix.  Any  system  of  the  form  of 
eq.  2.1  can,  by  proper  state  transformation,  be  put  into  the  form  of 
eqs.  2.2  and  2.3.  In  an  improper  system,  where  £2^  is  rank  deficient, 
the  output  measurements  which  track  the  command  inputs  are 


w  (t)  = 


C£l  +  M  ^1^  £2  +  M  (t)! 

j^X2  (t)j 

Cll'  F2]  fxi  (t)' 

X2  (t) 


(2.4) 


(2.5) 


where  M  is  an  extra  plant  measurement  matrix  which  must  be  designed  so 
that  £2  B2  has  full  rank.  The  trackers  examined  in  this  report  are 
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■'-iT  J  '-I  'JT  1 


I? 

b' 

I 

p: 

I 


improper,  requiring  a  measurement  matrix,  M.  The  state  vector  x(t)  s 
T 

(xj^  (t) ,  X2  (t))  in  these  equations  represent  the  kinematic  varicUsles 
used  in  the  aircraft  equations  of  motion 

Xi  (t)  -  ♦(t)  (2.6) 

X2  (t)  «  {0(t),  p(t),  r(t))’’  (2.7) 


i 

. 

!  C* 

•  •, 

i 


The  high-gain  error  actuated  analogue  controller  is  governed  on  T 
by  a  proportional  plus  integral  control  law  equation  of  the  form 

u  (t)  »  g(KQe(t)  +  l^(t))  (2.8) 

where  and  1^  are  controller  gain  matrices,  £  (t)  is  the  error  vector 

e  (t)  «  V  (t)  -  w  (t)  (2.9) 

V  (t)  is  the  command  vector,  eund  z  (t)  is  the  time  integral  of  e  (t) . 

2  (t)  =  e  (t)  -[v  (t)  -  w  (t)]  =[v  (t)  -  FtXi  (t)  -  F?xp(t)]  (2.10) 

A  functional  block  diagram  of  the  complete  tracker  is  shown  in  Figure 

2.1. 

It  is  desired  that  the  control  input  vector,  u(t) ,  cause  the  output 
vector,  ^(t) ,  to  track  any  constant  comm6uid  input  vector,  v^(t) ,  so  that 


lim  (V  (t)  -  ^  (t))  »  0  (2.11) 

t-*  <*> 

Such  tracking  is  possible  using  the  controller  form  of  eq.  2.8,  due  to 
the  fact  that  the  error  vector,  e  (t) ,  assumes  the  steady-state  value 
(Ref  7). 
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mii 


9 

t:-' 


i  CF 


n 


lim  e  (t)  =  lim  {v  (t)  -  w  (t) )  =0 
t  “*  <*  t  ® 


(2.12) 


The  steady-state  value  of  w  (t)  is 


lim  w  (t)  =  lim  (Fj^  +  £2  X2(t)) 


t  — *  ®  t 


(2.13) 


lim  ((Cl  +  M  ^1)3^  (t)  +  (C2  +  M  (2.14) 

t  -*  “> 

lim  ((£1  Xi(t)  +  £2  M(Aj^iXt  (t)  +  Ai  ^x-?  (t) ) ) 

t-»»  (2.15) 


=  lim  ( (CjL  XI  (t)  +  £2  2C2  (t) )  +  M  ^  (t) ) 


(2.16) 


t  -♦  » 


For  a  stable  system  with  a  constant  vector  input,  the  states  reach  con¬ 
stant  values,  and  therefore  it  is  evident  from  the  state  equation,  eq. 
2.2,  that 


lim  ^(t)  =  lim  (^ixx(t)  +^2X2(t))  =0 


(2.17) 


t— ♦  “  t  — ♦“> 


Therefore  eq.  2.16  yields 


lim  w(t)  =  lim  (£1  2^(t)  +  £2  X2(t)) 
t  — *  »  t  — ^  * 

=  lim  (Y.(t)) 
t  -»« 


(2.18) 


(2.19) 
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Thus,  in  the  steady-state  the  output,  ^(t)  tracks  the  command  input, 

.  v(t) ,  for  arbitrary  initial  conditions.  Equations  2.12  through^;  2.19 

« 

show  that  tracking  is  possible  for  any  measurement  matrix,  M,  provided 


that  the  condition  that  £2  £2  have  full  rank  is  satisfied. 

n-J.  S.  i 

In  the  edx3ve  equations,  (t)  e  R  ,  3^(t)  e  R  ,  ii(t)  e  R  , 

^(t)  e  R  ,  w(t)  e  R  ,  £(t)  e  R  ,  £{t)  e  R  ,  v(t)  e  R  ,  e  R  '  ' 


£1  e  R 


S.X  (n-i) 


£2  e  R 


M  e  R 


ix(n-i) 


The  closed- loop  equations  cure  formed  by  combining  eqs.  2.2,  2.3, 
2.5,  2.8,  and  2.10,  to  yield 


^(t)  =  0 

Lg 


g®2  ISl'  A?T-gB?  Kq  F]^,  A??-gB2  Kp  F2J 


2  (t) 


(t) 


=1^(0,  £1,  C2)J  rz  (t) 

2Cl(t) 

_3^(t) 


(2.20) 


(2.21) 


The  closed-loop  transfer  function  matrix  relating  ^■he  plant  output  vector 
to  the  command  input  vector  of  the  system  governed  by  eqs  2.20  and  2.21  is 


G(X)  =  (0,  Cl,  C2 


Ki 


-A^jl+gB^Kp  Fi, 


-^2 

-Xl£-^2+gB2KoF2. 


which  is  equivalent  to 


-1 


G(X)  =  CcL^XIn+£-  J  Bd, 


vAiere 


[:,i 

^  [  +®2  JSl'  ■  §2  lk>  l!l 

^  “  r  b22^  “  l2  Kp  £2] 


(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 


[ij 


By  letting  the  perturbation  parameter  e  =  1/g  go  to  zero  the  transfer 
function  matrix  ^(X),  approaches  (Ref  7)  the  asymptotic  form 


where 


auid 


r(X)  =  r(x)  +  r(X) 

r^(X)  =  ^  ^^in  “  ®D 
h\)  =  C2  [xi^  -  9  hAj  g  22 


(2.35) 

(2.36) 

(2.37) 


The  matrices  in  eq.  2.36  are 


r  -1  *1 

r  -1 

- 

[^^■22  A4  ^J  = 

-  Ko  Kl 

0 

-1 

-1 

A12F2  ^ 

A11-A12  22  2i_ 

(2.38) 


h2 


0 


h3 


]-D 


22^ 


1^,  Cl  -  £2  ^ 


(2.39) 


(2.40) 


The  "slow"  modes  of  the  tracking  system  are  defined  as  those  in  which 
the  eigenvalues  are  not  direct  functions  of  g.  These  are  the  poles  of 
(X)  and  are  easily  found  from 


r« 


Since  ^  is  lower  block  triangular,  eq.  2.38  yields  two  sets  of  slow 
poles  described  by 

z,  XeCtlX  K  +  K-  1  =0|  i  (2.42) 

and  4.  «  -o  j. 

*2  “  {  I  ‘  ^11  ^12  ^2*^  ^1  '  *  °  i 

Set  contains  the  transmission  zeros  of  the  system. 

The  "fast"  modes  of  the  tracking  system  are  those  in  which  the 
eigenvalues  eure  directly  dependent  on  g.  These  eigenvalues  eu:e  the 
poles  of  £(X)  and  are  easily  found  from 

det  [Xlj^  -  9  ^  ]  =0  (2.44) 

This  equation  yields  the  set  of  fast  poles  described  by 

Z3  =.jXeC:  |Xl^  +  9  £2  2-2  ^  I  “ 

In  order  for  tracking  to  occxur,  the  closed-loop  system  must 

obviously  be  stable,  requiring  that  all  the  eigenvalues  in  the  sets 

z^,  Z3  and  z^  be  in  the  left  half  of  the  complex  plane.  Therefore, 

for  sufficiently  Icurge  gains,  the  controller  matrices,  K  amd  K, , 

~o  ~1 

auid  measurement  matrix,  M,  must  be  chosen  so  that  z^^e  C”,  Z2  e  C  , 
and  Z3  e  C  . 

To  achieve  a  non-interacting  system  it  is  desired  to  have  a 
diagonal  asymptotic  transfer  function  matrix  £{X) .  in  order  for  the 
fast  asymptotic  transfer  function,  r^(X) ,  to  be  diagonal  it  is  neces- 

E.2  So  ■  <£2  *  !i  ^.12'  52  £0  -  i  =*  diag  a2t  aj,  (2.46) 

where  Oj  e  ( j  *  i,  2,  . . . ,  J,)  .  By  selecting  appropriate  values 
for  ojf  the  elements  of  can  then  be  determined 
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5o  -  ‘£-2  £2>  ■  ^ 


^us  the  fast  eusyiqptotic  eigenvalues  are 


*3  “  g»  -<^2 

^  may  be  chosen  ais  a  positive  multiple  of 


(2.47 


(2.48 


(2.49 


Solving  eq.  2.40,  the  first  set  of  slow  aisyii^totic  eigenvalues  are 


Si  -  {-h} 


(2.50 


There  sure  two  methods,  presented  by  Ridgley,  Baada.,  and  D'Azzo 
(Ref  9) ,  for  selecting  a  measurement  matrix,  M,  which  diagonalizes  the 
slow  asynptotic  transfer  function,  £(X).  Ihe  first  method  is  used  when 


n  >  21  and  v^en 


det  B*  =  det 


£l^  A’  B  • 


c^  A'^^l  B’ 


(2.51 


vdiere  A'  =  A^^^f  B'  *  ^2'  — '  “  — i'  £4  £' » 

d^  »  min  (j;  c^"^  A^^  B'  0,  j  *  1,  2,  ..,  n-1-1)  or  d^  =  n-l-l  if 

c^"^  A^^  B'  *  0  for  all  j . 

The  measurement  matrixr  M,  is  chosen  to  satisfy 


t,  •  ♦a*,,)  - r 


(2.S2 


(2.53) 


rr 


-  (=1  +  W  3il>  -  -  [ 

v^ere  6  is  the  maximum  value  of  d. ,  ^  cure  suit2d>ly  chosen  diagonal  matrices , 
euid 


kIo 


Nr  £'  ~ 


(2.54) 


It  is  not  always  possible  to  find  meeisurement  matrix,  M,  which 
satisfies  both  eqs.  2.52  and  2.53.  In  this  case  it  is  not  possible  to 
diagonalize  the  slow  eisynptotlc  transfer  ftmction  matrix,  r^(X). 

The  second  method  is  used  vdien  n  <  22.  and  det  »  0.  The  assignable 
elements  of  *  (jj^  +  M  ^2^ '  those  containing  elements  of  M,  are  per¬ 
mitted  non-zero  values  only  if  B*  has  a  non-zero  element  in  the  corre¬ 
sponding  position.  Also  should  be  chosen  to  be  diagonal.  Aged-n 

it  is  not  always  possible  to  find  a  matrix  M  which  diagonalizes  the  slow 
asymptotic  trcuisfer  function,  r^(A). 

When  there  is  no  measurement  matrix,  M,  which  will  decouple  the 
system  outputs  and  thus  make  the  asyn^totic  tremsfer  function,  r^(X) , 
diagonal,  then  M  must  be  iteratively  selected,  aiming  at  a  closed-loop 
system  that  is  st^d^le  and  nearly  non- interact! ve .  'Riis  is  a  limitation 
of  these  two  methods. 

When  it  is  desired  to  include  first-order  actuators  into  the  model, 
two  methods  are  availeQjle.  One  is  to  add  an  actuator  transfer  function 
to  the  forweurd  loop  path,  as  shown  in  Figure  2.2.  The  other  is  to 
augment  the  open- loop  state  equations. 
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.  The  actuator  transfer  function,  in  the  haplace  domain,  for  first- 
order  actuators  is  !! 


(s)  «  A  1  si  +  A 
—A  — C  I  — C 


(2.55) 


where  A  is  a  diagonal  matrix  contetining  the  I  values  of  positive  actuator 
constants.  The  resiilting  vector  differential  equation  is 


i^(t)  =  (t)  +  ^  u  (t) 


(2.56) 


Adding  the  actuator  transfer  fimction,  G^,  to  the  forweurd  path  augments 
the  closed- lo<9  equations,  eqs.  2.20  and  2.21.  Combining  eqs.  2.5,  2.8, 
2.9,  and  2.56  yields 


^(t)  -  -  ^  u^(t)  +  ^  ^g(K^  e(t)  +  £(t))J  (2.57) 

=  -he  +  9^  [Ko^H.<t)-Fj^Xj^(t)-F2X2(t)+I^z(t)]  (2.58) 


The  closed- loop  state  vector  is  augmented  to  include  ^  as  a  state.  The 
resulting  closed-loop  equations  are 


x^(t) 

X2(t) 

u^(t) 

—A 


”Li  '  ■L2  '  0 

-11  '  -12  '  ° 

-21  '  -22  '  -2 


<3h^'  -^Mo^l'  -^^c^2' 


z  (t) 


Xj^(t) 

X2(t) 

—A 


gA  K 
^-c-o 


(2.59) 


^(t)  »  |o»  £j^f  c^,  oj  r  (t) 


x^(t) 

X2(t) 

u^(t) 

—A 


(2.60) 
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To  add  the  first-order  actuators  by  augmenting  the  open- loop  state 
equations,  the  open-loop  state  vector  is  esq^emded  to  incltide  u^.  Tlie 
resulting  open-loop  equations  are 

.-21 '  -22  J  L-2'^^U  .-2  '- 


where 


2.(t) 


Xj^(t) 

X2(t) 


ii'  =  Fiitt)! 

-2'  = 


(2.62) 


(2.63) 


(2.64) 


-11'  "  -11 

A, 


21  -22 


(2.65) 


^2’  = 


U] 

[-^c] 

[  -1'  -4 

[5] 


(2.66) 


(2.67) 


(2.68) 


(2.69) 


(2.70) 


(2.71) 
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Equations  2.62  and  2.63  are  substituted  for  eqs.  2.2  amd  2.3,  then  the 
theory  outlined  in  this  chapter  is  applied. 

t 

In  this  report  the  first-order  actuators  are  added  by  augmenting 
the  closed-lo<^  equations. 


III.  Discussion  of  Results 


Design  Procedure 

t 

As  a  result  of  this  thesis  it  is  concluded  that  Porter's  design 
method  is  easy  to  apply.  Satisfactory  controllers  are  found  for  a 
specific  maneuver  and  flight  condition  without  many  design  iterations. 
Also  a  designer's  "feel"  is  quickly  obtained,  m£Ucing  controller 
designs  easy  to  obtain.  A  feel  is  when  the  designer  knows  how  differ¬ 
ent  parcuneters  affect  the  output  response.  The  effect  of  the  various 
parameters  discovered  during  the  design  iterations  for  the  con¬ 
trollers  of  this  report  aure  now  discussed. 

The  measurement  matrix,  M,  is  selected  as  a  sparse  matrix  using 
the  second  method  discussed  in  Chapter  II  (Ref  9) .  The  reciprocal 
of  the  non-zero  element  of  the  M  matrix  becomes  the  transmission 
zero  of  the  system.  A  transmission  zero  of  -4  is  selected  for  all 
of  the  controller  designs  in  this  report.  This  value  is  recommended 
by  Lt  Ridgely  (Ref  10) . 

The  value  of  the  sigma  matrix,  does  affect  the  outputs  in  a 
predictable  way.  To  keep  the  outputs  decoupled,  and  the  fast  asymptotic 

A 

transfer  function,  £,  diagonal,  ^  is  selected  as  a  diagonal  matrix. 

Each  diagonal  element,  a^,  affects  the  corresponding  output  (the  first 
diagonal  element  affects  the  first  output,  etc.).  Increasing  the 
value  of  the  element  causes  the  output  to  track  the  commemd  input 
closer.  A  large  diagonal  sigma  element  is  desirable  when  the  spe¬ 
cific  output  response  is  commanded  to  zero.  An  adverse  effect  of  a 
large  diagonal  element  is  that  the  magnitude  of  the  control  surface 
responses  are  increased.  Decreasing  the  diagonal  element  causes  an 


increase  in  the  transient  response  of  the  output.  Varying  the  diagonal 
elements  in  the  sigma  matrix  is  the  main  avenue  used  to  achieve  desired 
output  auid  control  surface  time  responses  for  the  controller's  designs 
contained  in  this  report.  An  observation  made  during  the  design 
iterations  is  that  an  optimal  sigma  matrix  found  for  a  specific 
maneuver  emd  specific  flight  condition  is  very  close  to  the  optimal 
sigma  matrix  for  the  same  meuieuver  at  a  different  design  flight  con¬ 
dition.  This  observation  is  shovm  by  comparing  the  sigma  matrices 
for  the  three  yaw  pointing  designs,  Tadsles  3-2  through  3-4. 

The  gain  value,  g,  is  used  in  the  design  iterations  for  the 
controllers  in  this  report,  to  adjust  the  peak  value  and  settling  time 
of  the  response  to  a  step  input.  The  gain  can  be  adjusted 

so  the  output  response  has  no  noticeadsle  peak  2uid  an  acceptable  settling 
time.  For  all  values  of  gain  an  overshoot  is  found  to  exist.  Very 
high  gain  and  very  low  gain  cause  a  large  overshoot.  The  gain  for 
a  minimum  overshoot  is  a  medium  value.  Settling  time  decreases  with 
increasing  gain  value. 

The  value  of  the  input  commands  Ccui  be  used  to  adjust  the  control 
surface  time  responses.  Changes  to  the  values  of  the  input  commands 
do  not  adversely  affect  the  output  responses,  in  that  they  will  track 
the  commands.  The  control  surface  responses  are,  naturally,  sensitive 
to  the  magnitude  of  the  input  commeuxd.  The  value  of  the  input  coramemd 
can  be  reduced  to  keep  the  control  surface  deflections  within  limits. 
This  is  the  last  parameter  to  be  adjusted  in  the  design  procedure. 

The  initial  values  of  input  ccxonands  eure  obtained  from  AFTI  wind 
tunnel  data  (Ref  12) . 
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The  control  surface  deflection  limits  for  the  designs  of  this  ' 


report  eure; 


Rudder,  -35°  5  5r  1  35° 

Aileron,  -35°  ^  i 
Cannard,  -30°  ^  Sc  1  30° 

These  limit  values  eure  selected  by  the  designer  and  considered 
reasoncdsle. 

The  interactive  computer  program,  listed  in  Appendix  C,  make 
rapid  design  iterations  possible. 


Robustness  Tests 

As  stated  in  the  Introduction,  the  robustness  tests  are  per- 
foamed  with  the  yaw  pointing  maneuver.  The  output  cat.rix  is 


(3.1) 


The  values  of  the  input  commands  to  produce  a  yaw  pointing  emgle  are 
shown  in  Ted)le  3-1.  The  side-slip  angle,  3,  is  commeuided  as  a  ramp  with 
the  slope  -3o  until  the  value  of  the  side-slip  angle  reaches  minus 
three  degrees.  At  this  time,  t,  the  slope  is  chcmged  to  zero.  The 
ramp  input  is  used  as  a  way  of  commeuiding  the  time  rate  of  the  side¬ 
slip  angle,  3.  The  roll  angle,  is  commanded  to  zero.  The  yaw  rate, 
r,  is  commanded  as  a  pulse  with  a  constant  magnitude  equal  to  3q.  The 
value  of  3o  is  a  function  of  the  dynamic  pressure,  Q,  which  is  de¬ 
termined  by  the  flight  condition.  At  time  t  the  value  of  yaw  rate 


is  commanded  to  zero. 


TABLE  3-1 


I 


i 

s. 


!  t* 

» 


INPUT  COMMANDS 


1 )  Side  Slip  Ang le, 

^(t)=-/8ot  0<t<T 

=  -  3  t  >  r 


2)  Ro  I  i  Ang  le,  0 

0(t)  :  0  t  >  0 

r 

3)  Yaw  Rate,  r 

t  < 

r(t)  =  0<t^r 


Q 

CONDITION 

SLOPE, ^ 

TIME 

CONSTANT.^ 

T 

0 

T 

0 

(deg/^ec) 

(seconds) 

(deg-sec) 

(degrees) 

Low 

•0.80 

3.75 

5.625 

3.00 

Low-Mid 

1,00 

3.00 

4.500 

3.00 

Mid 

■3.00 

1.00 

1.500 

3.00 

High 

4.00 

0.75 

1.125 

3.00 
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The  results  of  the  yaw  pointing  robustness  test  for  varying  flight 
conditions  eure  summcirized  in  Tables  3-2  through  3-4  and  Figures  3-2a 
through  3-2f.*  In  all  cases,  except  one  (Design  lA,  Q  »  825)  the  con¬ 
troller  performs  the  desired  yaw  pointing  maneuver.  This  is  shown 
by  the  fact  that  the  steady-state  side-slip  emgle,  8ss> 
steady-state  pointing  angle,  ilisst  ate  equal  in  magnitude  and  opposite 
in  sign.  The  pointing  cuigle,  \p,  is  the  time  integral  of  the  yaw 
rate 


il/Ct) 


r  (t)  dt 


(3.2) 


t 

The  most  noticeable  difference,  when  the  controller  is  used  outside 
its  design  flight  condition,  is  the  response  of  the  yaw  rate,  r. 
Outside  the  design  flight  condition  the  yaw  rate  has  an  initial  over¬ 
shoot  and  a  few  oscillations.  This  does  not  affect  the  performance 
of  the  maneuver,  based  on  the  smooth  response  of  the  side-slip 
angle  euid  the  pointing  angle. 

Figures  3-3  through  3-5  show  the  output  emd  control  surface  time 
responses  of  all  three  controllers  performing  a  yaw  pointing  maneuver 
at  a  specific  flight  condition.  Figure  3-3  shows  the  responses  of  a 
low  dynamic  pressure  maneuver  (Q  =  109) .  Figure  3-4  shows  the  re¬ 
sponses  of  a  medium  dynamic  pressure  maneuver  (Q  =  443) .  Figure  3-5 
shows  the  responses  of  a  high  dynamic  pressure  maneuver  (Q  <■  825) . 

The  controller  designed  at  a  low  dynamic  pressure.  Design  lA, 
consistently  overcommanded  the  control  surfaces  when  performing  the 
yaw  pointing  maneuver  outside  its  design  flight  condition,  as  shown 
•Figure  3-1  shows  the  basis  for  the  data  in  Tables  3-2  through  3-**. 
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in  Figure  3-2d.  Controllers  designed  at  medium  dynamic  pressure. 
Design  3A,  only  overcommanded  the  control  surfaces  at  two  of  the  ten 
flight  conditions  simulated,  and  the  value  of  the  deflection  error  is 
not  great. 

Figures  3-2e  and  3-2f  show  the  r2mge  of  motion  of  the  velocity 
vector  pointing  auigle.  The  velocity  vector  is  the  sum  of  the  side¬ 
slip  angle  and  the  pointing  amgle.  Figure  3-2e  shows  the  range  of 
values  that  the  velocity  vector  pointing  emgle  moves.  In  the  worst 
case  (Design  lA,  Q  =  825)  the  total  horizontal  translation  is  not 
greater  tham  10  feet,  during  the  time  of  the  maneuver,  t.  In  this 
case  the  steady-state  velocity  vector  pointing  angle  is  0.13  degrees, 
resulting  in  a  steady-state  horizontal  velocity  of  2.2  feet/second. 

Figure  3-2c  shows  the  time  integral  of  the  roll  angle.  This 
(Quantity  is  used  as  a  judgement  of  performance,  emd  a  basis  for 
comparison.  In  the  worst  case  (Design  3A,  Q  =  70)  the  value  is  0.176 
degree-second.  This  is  equivalent  to  a  roll  angle  of  0.06  degrees 
for  three  seconds.  The  magnitudes  of  the  roll  angle  time  integral 
shows  that  in  all  cases  the  pleme  remained  essentially  wings  level. 

Based  on  all  performance  criteria  shown  in  Tedjles  3-2  through 
3-4,  and  Figures  3-2a  through  3-2f,  the  controller  designed  for  a 
medium  dynamic  pressure  flight  condition.  Design  2A,  is  the  more 
robust  design. 

Another  robustness  test  is  judging  the  performance  of  a  con¬ 
troller  designed  to  perform  a  specific  mcuieuver  at  a  specific  flight 
condition  performing  other  memeuvers  at  the  same  flight  condition. 
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TABLE  3-2  Design  lA  Figures  of  Merit 

DESIGN  U  LOW  Q  Q  -  109  h  -  20,000  ft  Maoh 


0  10  0 

M  -  0.25  r »  0  15  0 

-  0  0  0  30 


20  0  0 
0  20  0 
0  0  20 


29.1768  1.6411  -6.4843 

1^-  5.4674  -4.3414  -1.6377  -  So  «  “  29.5 

44.6171  0.5608  20.4908 


70 


%error  tp  tg  tg^ 


2.993  1.0508  31.350  0.10  0.40  0.45 


109  -2.999  3.012  0.8197  2.462  0.35  i  0.40  0.30 


CONTROL 
4  SURFACE 
COMMEN 


exceeds 

0.0414 


Within 

0.0242  linits 


163  -2.999  3.005  1.0114 


0.10  0.45  0.40  0.0162 


223  -2.999  2.992  1.5377  53.770  0.05  0.40  0.45  0.0183 


436  -2  .999  3.018  3.7458  24.853  0.10  0.45  0.35  0.0168 


443  -2  .999  3.037  3.8850  22.820  0.15  0.45  P.35  0.0147 


532 


3.007  4.2513  41.710  0.05  0.40  0.40  O.OI53 


652  -2.999  2.9915  4.5585  51.883  0.05  0.40  0.40  0.0184 


825  -2.999  3.132  4.2150  5.375  0.60  0.75  p.lO 


0.0159 


198  -2.999  2.993  6.3388  58.470  0.05  0.40  0.15  I  0.0171 


exceeds 
limits 
>  10 


within 

limits 


exceeds 
linits 
Sc  >  23.9 


exceeds 
limits 
So  >  43.1 


exceeas 
limits 
>  23.1® 


exceeds 
limits 
S 


exceeds 
limits 
<fr  >10.6® 


exceeds 

limits 


^  - - o—  — 


DESIGTH  2A  Mp  Q  Q  -  443  Ib/f^  h  -  5|000  Mnoh 


* 

0 

;;i  «  0.25 
0 

m 


10  0 
0  15  0 
0  0  30 


20  0  0 
0  20  0 
0  0  20 


11,6944  0.1981  -2.0615 

2.7339  -1.0522  -0.1395  ^  S  -  29.0 

14.8517  -0.2256  4.4020 

•  « 


Q  ^ss 


70  -3.000  3.0U 


109  -3.000  3.001 


163  -3.000  2.995 


223  -3.000  2.5821 


436  -3.000  3.053 


443  -3.000  3.043 


532'  -3.000 1 3.020 


652  -3.000  3.011 


3.0801  2.6700 


825  -3.000  2.993  6.0860  52.1500  0.05  0.40  0.40  O.0367 


ts 

ts 

0.50 

0.45 

0.45 

0.45 

0.45 

0.45 

0,45 

0,40 

0.30 

0.30 

0.20 

0.35 

0.40 

0.40 

0,50 

0.35 

0.40 

0.40  ^ 

vlthin 

liDitS 


v,t-thin 

liaita 


vithin 

liaita 


Vi  thin 
liaita 


vithin 

liaita 


vithin 

liaita 


vithin 

linita 


vithin 

liaita 


98  -3.000  2.981  5.4505  36,2625  0,05  0.45 


0.0365 


TABLE  3-^  Design  3A  Figures  of  Merit 

DCSIGH  3A  high  Q  Q  «  825  Vbft^  L  .  10,000  ft  Mach 


M»  0.25 
0 


*2  0  0 
£«  0  15  0 
0  a  30 


20  0  0 

4  =1  0  20  0 

0  0  20 


22.5871  0,1857  -1.6393 

5.9893  -0.7110  -0.0607  £.1“  %  e  =  29.5 

22.3193  -0.1024  1.9477 


Q  P 

^  ss 


70  -3.000 


109  -3,000 


3.045  1.1625  45.312 


0,80  ;o.6C  [  0.1763 


3.025  1.2601  57.512  0.05  |0.50  0.55  0.1119 


-3.000  3.021  1.3857 


0.05  0.50  0.45  0,0829 


223  -3.000  3.009  1.6871  j  68,710  j  0.05  0,45  0.50  0,1361 


436  -3.000  2.972  4.1703  39.010  0.10  0.50(0,60(0.0705 


443  -3.000  2.952  4.0944  36.480  0.10  0.40  0.65  0,0662 


532  -3.000  2.907  3.6494  21,646  0.15  0.95  0.60  0.0632 


552  -3,000  2.951  3.0646  2,153  0.30  0,40  jo, 55  0,0652 


825  -3.000  3.041  4.1046  2.615  0.35  0.45  0.60  0.0635 


198  -3.000  3.045  5.4566  36.415  0.10. 
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SURFACE 
COM  MEN 


within 

linits 


within 

lioits 


within 

lixaits 


within 

linits 


exceeds 
linits 
oc  >8.6® 


exceeds 

linits 


within 

limits 


within 

linits 


0.0635 

within 

1 ini ts 

.55 

0,0655 
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limits  * 

I  DEGREE- SECONDS 
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OUTPUTS  OUTPUTS  OUTPUTS 


OUTPUTS  OUTPUTS  OUTPUTS 


OUTPUTS  OUTPUTS  OUTPUTS 


The  controller  design  2A  is  given  ccnuncuids  to  perform  a  wings  level 
turn,  a  horizontal  translation,  and  a  constant  altitude  coordinated 
turn.  The  magnitudes  of  these  ccmtmands  are  adjusted  until  the  maneuver 
can  be  performed  without  commanding  the  control  surfaces  beyond  their 
limits . 

The  input  command  for  a  wings  level  turn  is  a  step  input  for  yaw 
rate,  r,  while  side-slip  angle,  3,  and  roll  angle,  (p,  are  commcUided 
to  zero.  For  horizontal  translation  the  side-slip  eUigle  is  commanded 
as  a  ramp,  till  a  side-slip  cUigle  of  four  degrees  is  reached.  The 
roll  cUigle  emd  yaw  rate  are  commanded  to  zero.  For  a  coordinated 
turn  the  relation  between  yaw  rate  and  roll  angle  is 

r  .  pi”  ♦  (3.3) 

while  side-slip  angle  is  commanded  to  zero. 

The  output  and  control  surface  time  responses  of  Design  2A 
perfoming  the  above  three  maneuvers  are  shown  in  Figure  3-6.  These 
responses  ^u:e  ccmipcired  with  those  of  Figure  3-7  which  show  the 
responses  of  controllers  designed  to  specifically  perform  the  maneuvers 
at  the  design  2 A  flight  condition  (Q  =  443) .  It  is  easily  seen  that 
Design  2A  can  not  perform  the  maneuvers  as  well  as  the  design  for  that 
specific  maneuver.  The  magnitudes  of  the  maneuvers  of  Design  2 A  are 
much  less  thcui  the  magnitudes  of  the  Mid  Q  designs.  For  example  the 
Mid  Q  design  for  horizontal  treuislation  can  successfully  commeuid  a 
side-slip  angle  of  four  degrees  in  1.33  seconds,  while  the  controller 
of  Design  2A  took  4.5  seconds.  In  this  respect,  the  controller  is 
not  robust  with  respect  to  varying  commemd  inputs.  The  Mid  Q 
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Figure  3-T  Kid  Q  Additional  Maneuvers 


Table  3-5  Mid  Q  designs 


MID  Q  WINGS  LEVEL  TURN 
M  =  r  =  fo 


.500 
0  10  0 
0  0  .5 


[So  0  0 

0  20  0 
0  0  20 


r  5.8472  0.1321  -0.03441 

IC=  I  1.3695  -0.7015  -0.00231  Ki  =  ^  g  =  20 
I  7.3259  -0.1504  0.07341 


Desgign  2A 
Mid  Q 


Steady-state 
Yaw  Rate 

Peak 

Overshoot 

1.20  deg/ sec 

1.75  deg/ sec 

2.0  deg/ sec 
1.8  deg/ sec 

MID  Q  HORIZONTAL  TRANSLATION 

[0  1  r  5  0  o' 

0.251  X  =  I  0  50  0 

0  ]  I  0  0  75 


?  50  0  4=  0  20 


20  0  0 
0  20  0 
0  0  20 


58.4719  0.6604  - 

13.6946  -3.5073  - 
73.2587  -0.7521  1 


-5.1539 

-0.3487  = 

11.0051J  ^ 


Design  2 A 
Mid  Q 


MID  Q  COORDINATED  TURN 


Steady-state 
Side-slip  angle 

Maneuver 

Time 

4.00  degrees 
4.00  degrees 

4.50  seconds 

1.33  seconds 

r  0  1  f  1  0  0] 

M  =  I0.25  r  =  0  1  0  A  = 

-  I  0  J  I  0  0  j 

(11.6944  0.0132  -0.0697] 

2.7389  -0.0701  -0.00461  K  = 
14.6517  -0.0150  0.14571 


l^=  2.7389  -0.0701 

14.6517  -0.0150 


Steady-state 
Yaw  Rate 


Design  2A 

6.15  degrees 

Mid  Q 

20.93  degrees 

r  20 

0 

d 

A  =1  0 

20 

0 

-  1  0 

0 

20 

g  = 

20 

controller  designs  and  performance  comparison  summary  are  shown  in 
Table  3-5. 
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Actuator  Dynamics 

As  discussed  in  Chapter  2  there  are  two  methods  to  add  first-order 
actuator  dynamics  to  the  system.  The  method  used  in  this  report  is  to 
augment  the  closed- loop  system  matrix.  This  does  not  cause  the  size 
of  the  measurement  matrix  to  increase.  If  the  actuator  dynamics  are 
added  by  augmenting  the  open-loop  system  matrix,  the  size  of  the 
measurement  matrix  increases  from  a  3-by-l  matrix  to  a  3-by-4  matrix. 
When  an  attempt  is  made  to  add  the  actuator  dynamics  this  way,  the 
method  used  for  selecting  a  measurement  matrix  that  both  stabilizes 
the  system  and  diagonalizes  the  asymptotic  transfer  function  fails. 

This  leaves  the  selection  of  a  measurement  matrix  to  a  trial  and 
error  iterative  procedure. 

The  actuator  dynamics  are  significant  in  the  simulations.  For 
comparison  the  yaw  pointing  robustness  test  is  repeated  for  the  same 
three  controllers  without  actuator  dynamics  in  the  model.  The 
results  of  these  simulation  tests  are  very  misleading.  Figures  3-8 
through  3-10  show  the  output  and  control  surface  time  responses  of 
the  three  controller  designs  performing  the  yaw  pointing  maneuver  at 
a  specific  flight  condition.  When  these  figures  are  compared  with 
Figures  3-3  through  3-5  the  significance  of  the  actuator  dynamics 
is  clearly  shown.  Note  that  in  a.  1  cases  the  outputs  follow  the 
cornmcmds  exactly,  with  no  overshoots  or  oscillations.  Also  note  the 
Smooth  control  surface  response  curves.  The  actuator  dynamics  in  the 
model  greatly  affect  the  time  responses  euid  should  not  be  ignored 
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when  designing  a  controller.  This  same  trend  is  observed  when 
Figure  3-11  is  compared  with  Figure  3-6. 

Without  actuator  dynamics  the  output  responses  are  much  jnore 
insensitive  to  changes  in  the  diagonal  elements  in  the  sigma  matrix  and 
gain  values.  This  is  another  reason  why  the  actuator  dynamics  should 
not  be  ignored  when  designing  a  controller. 

Also  studied  is  the  effect  of  changing  the  first-order  actuator 
dynamics  time  constant.  The  controller  is  designed  with  all  actuator 
constants  at  20.  This  design  is  used  in  the  simulation  of  a  yaw 
pointing  maneuver  with  the  actuator  constants  varying.  The  responses 
are  shown  in  Figures  3-12  through  3-18.  The  three  numbers  on  the 
title  of  the  response  correspond  to  the  three  actuator  constauits. 

The  respective  order  is;  rudder  actuator  constemt,  aileron 
actuator  constant,  and  cannard  actuator  constant.  For  example, 

(20,  10,  40) represents  a  rudder  constemt  of  20,  a  aileron  constant 
of  10,  and  a  camnard  constant  of  40. 

Fig\ire  3-12  shows  the  responses  when  all  constants  are  de¬ 
creased  and  increased  by  a  factor  of  two.  Figure  3-13  shows  the 
responses  when  just  one  constant  is  decreased  by  a  factor  of  two. 

Figure  3-14  shows  the  responses  when  just  one  constant  is  increased 
by  a  factor  of  two.  Figure  3-15  shows  the  responses  when  two  constamts 
eure  decreased  by  a  factor  of  two.  Figure  3-16  shows  the  responses 
when  two  constamts  are  increased  by  a  factor  of  two.  Figure  3-A7 
shows  the  responses  when  one  constant  is  decreased  and  one  constamt 
is  increased  by  a  factor  of  two. 


This  actuator  study  shows: 


1.  The  aileron  has  little  effect  in  the  yaw  pointing  maneuver 

2.  Decreasing  the  cannard  constant  increases  settling  time  and 
increases  cannard  activity 

3.  Increasing  the  rudder  constant  causes  increased  rudder 
activity 

4.  Increasing  both  rudder  and  cannard  constants  cause  the 
magnitude  of  surface  deflection  to  decrease. 

Overall  this  study  amplifies  the  importance  of  actuator  dynamics  in 
the  model,  in  that  they  significantly  affect  the  output  and  control 
surface  time  responses.  The  actuator  time  constant  is  dependent  on 
the  physical  system  used  in  the  aircraft.  Varying  the  constants  in 
the  design  can  be  advantageous  in  a  preliminary  design  study  and  can 
be  used  as  a  basis  for  selecting  the  best  actuator. 

In  this  study  there  is  no  constraint  placed  on  the  rate  of  con¬ 
trol  surface  deflection.  As  a  result  the  control  surface  deflection 
rates  were  often  unreasonably  fast.  Any  further  study  should  include 


this  constraint. 
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IV.  Conclusions  and  Recommendations 


Conclusions 

Major  conclusions  are  based  on  the  results  of  this  report.  They 
are  first  presented  in  Chapter  Ill  and  are  now  summarized. 

The  design  method  of  Dr  Porter  is  easy  to  apply  and  does 
produce  controller  designs  to  perform  desired  maneuvers.  This  method 
is  relatively  new  and  promising.  Research  and  applications  of  this 
method  should  continue. 

The  controllers  eire  robust  with  respect  to  varying  flight  con¬ 
ditions.  They  are  not  robust  with  respect  to  varying  maneuver  com¬ 
mands.  The  controller  designed  at  a  medium  dynamic  pressure  flight 
condition  produced  the  more  robust  controller  design. 

Actuator  dynamics  significantly  affect  the  output  and  control 
surface  time  responses.  Simplified  models  without  actuator  dynamics 
included  provide  misleading  results.  Actuator  dynar,'.ics  should  not 
be  ignored  in  the  design  of  flight  controllers  with  this  method. 

Recommendations 

This  study  followed  the  study  performed  by  Lt  Ridgely  (Ref  10) . 
His  study  examines  the  longitudinal  modes  of  AFTI .  This  study 
examines  the  lateral-directional  modes  of  AFTI.  The  next  study  should 
combine  the  two  modes  and  examine  the  response  of  the  total  aircraft 
when  performing  various  maneuvers. 

Further  robustness  examination  should  be  conducted  to  find  if  a 
controller  cem  be  designed  that  is  both  robust  with  respect  to  flight 
condition  euid  maneuver  commands.  If  so,  what  parameters  affect  ro¬ 
bustness. 
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The  computer  program  used  for  this  study  is  limited  when  compared 
with  the  computer  program  developed  for  the  design  of  digital  con¬ 
trollers  using  Porter's  method,  MULTI  (Ref  13) .  MULTI  should  be 
augmented  to  include  the  design  of  emalogue  controllers. 

The  precision  of  program  ZERO  should  be  increased  from  single 
precision  to  double  precision  to  increase  the  program's  accuracy. 

Also  the  subroutines  in  ZERO  that  check  the  rank  of  the  system 
matrices  should  be  changed  so  that  complex  system  zeros  can  also  be 
tested. 
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Appendix  A-1.  Development  of  the 
Aircraft  Equations  of  Motion 


As  stated  in  the  introduction,  only  the  lateral-directional,  de¬ 
coupled  maneuvers  eu:e  analyzed  in  this  report.  Hence,  only  the  lateral- 
directioneU.  set  of  equations  are  developed  here.  The  airfreuae  model  is 
a  convention2d. ,  small  perturbation  set  of  equations,  which  are  lineeurized 
ad>out  straight  and  level  flight  (6^  »  .  All  four  m£uieuvers  emalyzed 

are  described  by  body  axis  varied^les,  therefore,  the  equations  of  motion 
eure  body  £ucis  equations.  The  linearized  equations  in  the  Laplace  domain 
are  (Ref  5) : 


sv  +  U^r  -  -  g  cos 

=  y  V  +  Y  p  +  y  r  +  ZY.d  (A.l) 

V  p  r  6 

sp  =  L'  V  +  L*  p  +  L*  r  +  EL'  6  (A.2) 

V  p  r  0 

sr  =  N*  V  +  N*  p  +  N*  r  +  EN'  ,6  (A. 3) 

V  p  r  6 

S(^  3  p  +  r  teui  6  (A.4) 

o 


These  equations  are  tremsformed  to  the  time  domain,  emd  arremged  to  have 
the  matrix  form 


X  =  A  X  +  B  u 


where 


T 


X  » 


T 


u 


(A.  5) 


(A.6) 

(A.7) 
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0 

,  0. 

1 

t 

taui  6 

o 

A 

a 

g  cos  6 

A 

'  ^v' 

(Y  +  w 
P 

o)/u. 

»r  - 

(A.  8) 

9 

,  L'p, 

f 

‘■’r 

Q 

N* 

P 

f 

N* 

r 

• 

0 

i 

0 

»  0 

B 

= 

^*6r 

f 

^*6a 

(A.  9) 

^’6r 

t 

^'6a 

**'6r 

t 

”'6a 

'  ***60 

The 

AFIT 

-16  data 

(Ref  12) 

is  not  consistent  with  eqs. 

A.l  through 

A.  9,  since  the  stability  coefficients  are  in  non-dimensional  form  amd 
are  for  the  stadaility  axis  system.  Also  some  of  these  stability  coef¬ 
ficients  have  units  of  per  radiaui,  while  others  have  units  of  per  degree. 
The  moments  of  inertia  given  are  in  dimensional,  body  axis  form. 

The  steps  used  to  transform  the  given  data  to  data  consistent  with 
eqs.  A.l  through  A.  9  aure: 

1.  Trauisform  all  stability  coefficients  with  units  of  per  degree 
to  per  radiam. 

2.  Trauisform  the  moments  of  inertias  to  the  stability  axis 
system. 

3.  Diraensionalize  the  stability  coefficients. 

4.  Trauisform  the  L  and  N  dimensional  stability  coefficients  to 
prime  form. 


6] 


5.  Treuisform  the  dimensional  stability  coefficients  to  the  body 


eocis  sytem. 

Ten  different  flight  conditions  are  used  in  this  study,  reqviirlng 
ten  sets  of  equations  of  motion,  ^e  data  for  these  flight  conditions 
are  listed  in  Tzdsle  A-1.  A  coi^uter  program  is  written  to  perform  all 
of  the  cibove  transformations.  Use  of  the  computer  program  guareuitees 
that  the  data  transformations  eure  applied  consistently.  An  exeui^le  euid 
program  listing  2u:e  in  Appendix  A-2. 

The  stability  coefficients,  as  originally  given,  in  non-dimensional 
st2d>ility  cucis  form  £u:e  listed  in  Table  A-2.  The  dimensionail  tremsforma' 
tlon  of  step  1  consists  of  multiplication  by  180°/ir.  The  trims  formation 
of  step  2  makes  the  multiplication  factors  of  steps  3  and  4  consistent, 
that  is,  inertias  emd  ctability  coefficients  are  all  in  the  stid>llity 
axis  system.  The  transformation  from  body  axis  to  sted>ility  axis  is 
accomplished  by  using  the  following  equations: 


(I  ) 

(I  K 

2 

cos 

9  - 

2(1  ). 

cos  e 

sin  e  ^ 

XX  s 

XX  b 

o 

xz  b 

o 

o 

+ 

t^zz^b 

,  2 
sin 

0 

O 

(A.  10] 

(I  ) 

cos^ 

2(1  ). 

cos  e  ^ 

sin  0 

zz  s 

zz  b 

o 

xz  b 

o 

°  o 

+ 

"xx'b 

.  2 
sin 

®o 

(A. 11] 

'  "Wb  -  “zz’b>  9o  ^ 

+  (I  ).  (cos^e  -  sin^e  )  (A.  12) 


TABLE  A-1 


Flight  Condition  Data 


Dynamic 
Pressure,  Q 

Altitude 

Density, p 

Sonic 

Velocity 

Mach 

Velocity ,  U 

Angle  of 
Attack , 

(lb/ft2) 

(feet) 

(slug/ft^) 

(ft/sec) 

(ft/sec) 

(degree ) 

70 

30,000 

0.000889 

994.6 

0.4 

397.84 

11.400 

109 

20,000 

0.001266 

1036.8 

0.4 

414772 

6.952 

163 

10,000 

0.001755 

1077.4 

0.4 

430.96 

4.552 

223 

40,000 

0.000585 

968.1 

0.9 

871.29 

3.692 

436 

20,000 

0.001266 

1036.8 

0.8 

829.44 

2.285 

443 

5,000 

0.002048 

1097.1 

0.6 

658.26 

1.824 

532 

30 

0.002377 

1116.4 

0.6 

669.84 

1.547 

652 

10,000 

0.001755 

1077.4 

0.8 

861.92 

1.733 

825 

10,000 

0.001755 

1077.4 

0.9 

969.66 

1.462 

1198 

30 

0.002377 

1116.4 

0.9 

1004.76 

1.247 

The  multiplication  factors  required  to  dimensionalize  the  stediility  coef¬ 
ficients,  step  3,  are  shown  in  Tcdsle  A-3.  l^e  trws formation  of  step  4 
is  acccm^lished  by  using  the  equations 


,  L.  +  (d  )  /(I  )  )N. 

'  =  I  xz  s'  XX  X  1 


(A. 13) 


(1  -  ((I  )„/(!  ),(I  ) J 

XZ  S  XX  s  zz  s 


'i  -  <«xz>^<V.«zz'.> 


(A.  14) 


The  fixed  adrcraft  parameters  required  for  steps  3  auid  4  for  the  AFIT-16 
are  shown  in  Table  A-4.  The  transformation  from  the  stjdsility  £ucis  system 
to  the  body  eucis  system,  step  5,  is  accomplished  with  the  following  set 
of  equations: 


(A.  15) 


<Vb  “  ^^p^  ‘  ^Vs  ®o 


(A. 16) 


<Vb  =  <^r)s  ®o  +  <^p>s 


sin  0, 


(A. 17) 


(L’  x).  =  (I*'  -)  cos  0  -  (N*  ,)  sin  0 

v,o  b  v,o  s  o  v,o  s  o 


(A. 18) 


<L'p)b  =  <^'p>s  ®o  ’  <^^'r>s  +  ®o  ®o 


+  (N'p)^  sin  0^ 


tb'r)b  =  (b'Pg  cos  0Q  -  ((N*j.)g  -  (L*  )j,  sin  8^,  cos  8^ 


-  (N*  )  sin2  0 
p  s  o 


(A.  19) 


(A. 20) 


<**'v,6’b  =  <N'v,6>s  ®o  +  <^’v,6>s  ®o 


(A.21) 
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=  (N-p)^  cos" 

0  - 
o 

((N*  )„  -  (t*  ),) 
r  s  p  s 

sin  0  cos  0 
o  o 

- 

(L*  )  sin^  9 
r  s  o 

(A.22) 

=  (N'^)^  cos^ 

6  + 
o 

((L'r),  +  (N-p)^) 

sin  0^  cos  0^ 
o  o 

+ 

(J.'p)g  sin^  0^ 

(A.23) 

The  side-slip  euigle,  g,  is  defined  as 

S  =  v/u  (A.  24) 

The  dimensioned.  3  sted>ility  coefficients  differ  from  the  v  stability  coef' 

ficients  by  a  factor  of  l/fj.  Note,  for  exeunple,  that  N.  3  =  N  v.  The 

p  r 

y*  terms  represent  a  multiplication  by  l/u»  that  is 

Y*i  =  (A.  25) 


Yhe  matrices  A  and  B  for  ed.1  ten  flight  conditions  are  listed  in  Tedsle 


^  TABLE  A-2 

AFTI-16  Non-Dimensional  Stability  Axis  Data 
Lateral-Directional  CN 


Q 

c 

C 

C 

C 

C 

C 

"6 

n 

P 

n 

r 

"Sdt 

2 

-1 

-1 

-1 

-1 

-1 

-1 

(Ib/ft 

) (deg  ) 

(rad  ) 

(rad  ) 

(deg 

(deg 

(deg 

70 

.002681 

-.038417 

-.496924 

-.001667 

.000379 

-.000979 

109 

.002121 

-.020725 

-.492522 

-.001618 

.000141 

-.000986 

163 

.001784 

-.013125 

-.487107 

-.001557 

.000007 

-.000955 

223 

.001800 

-.010649 

-.279610 

-.001436 

-.000084 

-.001069 

436 

.001710 

-.006850 

-.346045 

-.001333 

-.000144 

-.000978 

443 

.001775 

-.005461 

-.476796 

-.001379 

-.000137 

-.000893 

532 

.001761 

-.004666 

-.473181 

-.001342 

-.000147 

-.000862 

^52 

.001677 

-.005165 

-.341586 

-.001229 

-.000149 

-.000899 

825 

.001636 

-.004298 

-.273077 

-.001137 

-.000149 

-.000885 

1198 

.001447 

-.003620 

-.260486 

-.000947 

-.000121 

-.000741 

Lateral-Directional 

,  CL 

Q 

C, 

C, 

C, 

C, 

C, 

C, 

1 

P 

1 

r 

^6r 

^6dt 

2 

-1 

-1 

-1 

-1 

-1 

-1 

)  (deg  ■^) 

(rad  ■*■) 

(rad  ) 

(deg  ) 

(deg 

(deg  ■'■) 

70 

-.002583  -.223766 

.164036 

.000179  -.001508 

-.001857 

109 

-.002396  -.240803 

.103189 

.000323  -.001899 

-.001649 

163 

-.002162  -.242331 

.064776 

.000387  -.002032 

-.001673 

(deg 

.000934 

.001071 

.001118 

.001196 

.001215 

.001168 

.001169 

.001217 

.001250 

.001228 


(deg 


.000014 

.000072 

.000127 


223  -.002358  -.362413  .087114  .000331  -.002195  -.001836  .000226 


TABLE  A-2  (cont) 


Q 

^1 

P 

^1 

r 

Cl 

6dt 

{lb/ft‘ 

(deg’^) 

(rad 

(rad 

(deg 

(deg 

(deg 

436 

-.001929 

-.313529 

.031455 

.000347 

-.001971 

-.001810 

443 

-.001703 

-.232873 

.003219 

.000392 

-.002030 

-.001699 

532 

-.001608 

-.230739 

-.004990 

.000386 

-.001954 

-.001656 

652 

-.001738 

-.306430 

.011514 

.000326 

-.001738 

-.001702 

825 

-.001742 

-.342130 

.006361 

.000294 

-.001523 

-.001682 

1198 

-.001602 

-.328623 

-.002132 

.000237 

-.001184 

-.001542 

Lateral-Directional  CY 

Q 

C 

C 

C 

C 

C 

C 

^r 

^6dt 

(lb/ft‘ 

(deg‘^) 

(rad 

(rad 

(deg 

(deg“^) 

(deg 

70 

-.002015 

.224381 

.603623 

.003064 

.000703 

.002916 

109 

-.021784 

.189589 

.557355 

.003069 

.000667 

.002528 

163 

-.020995 

.141151 

.535556 

.003066 

.000571 

.002132 

223 

-.022820 

.110272 

.556043 

.002667 

-.000138 

.001998 

436 

-.021901 

.070000 

.547367 

.002622 

-.000113 

.001618 

443 

-.021095 

.058391 

.537543 

.002991 

-.000080 

.001364 

532 

-.020953 

.051413 

.537411 

.002920 

-.000094 

.001280 

652 

-.021466 

.056100 

.543785 

.002422 

-.000112 

.001463 

825 

-.021582 

.045047 

.548795 

.002110 

-.000077 

.001492 

(deg"^) 

.000205 

.001409 

.000144 

.000199 

.000224 

.000217 

C 

^60 

(deg'^  ) 

.001006 

.001075 

.001233 

.001350 

.001509 

.000146 

.001462 

.001622 

.001785 


1198  -.020906  .038223  .537804  .001745  -.000063  .001208  .001793 


TABLE  A-3 


Lateral-Directional  Dimensionalization  Factors 


2m  y3 

II 

CO. 

pSU^b 

2Izz 

n 

4m  yr 

N  = 
r 

pSUb^ 

4Izz 

^nr 

L 

r 

3b  p 

4m  yP 

Np  = 

pSUb2 

^^zz 

f^np 

^  c 

2m  y« 

Nfi  = 

pSU^b 

2Izz 

Lj  =  £ 

TABLE  A-4 

AFTI-16  Fixed  Aircraft  Parameters 


weight,  w 

21018  pounds 

spam,  b 

30  feet 

Area,  S 

300  (feet)^ 

(Ixx)b 

10033.429  slug-ft 

(Izz)b 

61278.452  slug-ft 

282.132  slug-ft 

TABLE  A- 5 

APTI-16  Plamt  Matrices 

X  =  Ax  +  Bu,  X  ^  =  [((),  3,  p,  rj  u  =  1^6^,  6^/  6^,] 
Design  1  Q  =  109 


0.0000 

0.0000 

1.0000 

.1219 

0.0000 

0.0000 

0.0000 

.0770 

-.1506 

.1210 

-.9926 

.0212 

.0090 

.0074 

0.0000 

-14.7042 

-  .8989 

.4595  ®  " 

2.8507 

-32.7814 

-.2979 

0.0000 

1.5960 

.0007 

-.2750 

-1.4230 

-  .4112 

.9809 

Design  2  Q  =  443 


0.0000 

0.0000 

1.0000 

.0318 

.0489 

-.3744 

.0318 

-.9995 

0.0000 

-39.9264 

-2.1137 

.0797 

0.0000 

6.2376 

-.0063 

-.7074 

■  0.0000  0.0000  o.oooo' 

.0531  .0046  .0250 
9.7929  -55.7388  2.6033 
-5.0544  -1.8927  4.3880 


Design  3  Q  =  825 


0.0000 

0.0000 

1.0000 

.0255 

.0332 

-.4836 

.0255 

-.9997 

0.0000 

-75.3125 

-3,9178 

.0382 

0.0000 

10.6995 

-  .0294 

-.5117 

Q  =  70 

0.0000 

0.0000 

1.0000 

.2016 

.0793  -.1025  .1977  -.9803 
0.0000  -11.0376  -.6157  .4996 
0.0000  1.2029  .0009  -.1777 


'  0.0000  0.0000  0.0000 

.0473  .0066  .0399 

13.4743  -82.0763  8.3887 

-7.7774  -  3.2921  8.7542 


0.0000  0.0000  0.0000 

.0143  .0067  .0047 

1.7994  -7.0917  -.7020 

-.9382  -  .1855  .5372 
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TABLE  A- 5  (cont) 


Q  =  163 

0.0000  0.0000  1.0000  .0796 


0.0000 

-19.1782 

-1.2617 

0.0000 

2.1153 

-.0009 

Q  =  223 

0.0000 

0.0000 

1.0000 

.0369 

-.1532 

.0644 

0.0000 

-28.1034 

-1.2593 

0.0000 

2.9444 

-.0149 

Q  =  436 

0.0000 

0.0000 

1.0000 

.0388 

-.3029 

.0399 

0.0000 

-44.5171 

-2.2233 

0.0000 

5.7767 

-.0167 

Q  =  532 

0.0000 

0.0000 

1.0000 

.0480 

-.4392 

.0270 

.2752 

-.1532 


.039S 


0.0000  -45,1440  -2.4708 

0.0000  7.4966  -.0081  -.83 


0.0000 

5 


4.2081 

2.0671 


0.0000  0,0000 


-20.3236  .3607 

-.6771  1.5439 


0.0000 


9 


4.7528 

2.6169 


0.0000  0.0000 


-30.0084  1.7593 

-1.1128  2.2662 


0.0000 

63 


0.0000 


0.0000 

9 


8.8177 

-53.9258 

3.6280 

4.7193 

-2.0255 

4.4971 

0.0000 

0.0000 

0.0000 

.0612 

.0047 

.0306 

1.4009 

-64.6667 

3.2286 

5.9204 

-3.2106 

5.2760 

70 


TABLE  A-5  (cont) 


Q  =  652 


0.0000  0.0000  1.0000 

373  -.4275  .0302 

0.0000  -59.6596  -3.1213 


0.0000 


8.6325  -.0212  -. 


0.0000 

.0482 


0.0000 

.0051 


11.9760  -72.1512 


-6.6298 


-.7405 


0.0000 


5.6206 

6.7320 


Q  =  1198 


0.0000 


0.0000  1.0000 


-.6574 


.0218 


0.0000  -100.3136  -5.2784  -. 


0218  0.0000 
9998  g  ^  .0549 

0773  15.6158 


0.0000 


13.7921  -.0377 


0.0000 

.0075 


-9.4350  -3.880 


0 . 0000 
.0564 


-96.4546  12.0815 


12.49861 
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THIS  PROGRAM  CONVERTS  NON-DIMENSIONAL-STABILITY  AXIS- 
DATA  TO  DIMENSIONAL-BODY  AXIS-DATA  FOR  THE  AFTI-IB 

SELECT  MODE.  LONGITUDNAL  ONLY  (1) 

LATERAL  ONLY  (2) 

BOTH  (3)  :  2 

STANDARD  AFTI  AIRCRAFT? 

YES=1.  N0=0.  :  1 

FIXED  AIRCRAFT  PARAMETERS 

WEIGHT  (LBS)  :21018. 

SPAN  (FT)  :30. 

AREA  (FT  **2) :300. 

MAC  (FT)  Ill. 32 

I-X  110033.429 

I-Y  153876.269 

I-Z  161278.452 

I-XZ  1282.13217 


INPUT  OARING  PARMETERS 

HEIGHT  (FT)  I  5000 
MACH  I  0.6 

SONIC  VELOCITY  (FT/SEC)  I  1097.1 
DENSITY  (SLUG/FT  »-*3)  I  0.002048 
ANGLE  OF  ATTACK  (DEG)  I  1.824 

MOMENT  COEF. 

CL  :  .1572 
CM  :  ,0026 
CD  :  .0231 

INPUT  THE  LATERAL  STABILITY  COEFICIENTS 
CN  ,  CL  ,  CY 

BETA  (DEG)  I  . 00 1 775 , - . 00 1703 , - . 02 1 095 
P  (RAD)  :  -.005461 ,-.232873, .058391 
R  (RAD)  :  -.476796, .003219, .537543 
RUDDER  (DEG)  I  -. 00 1 379 ,. 000392 00299 1 
FLAPERON  (DEG)  I  -. 000 137 ,-. 002030 ,-. 000080 
DIFFAIL  (DEG)  I  -. 000893 ,-. 00 1699 ,. 00 1 364 
CANNARD  (DEG)  I  . 00 1 168 , . 000 1 46 , . 00 1 409 
INPUT  ALERON  COMBINATION  TERMS 
ALERON=  A^FLAPERON  +  B*DIFF  TAIL 
ENTER  A,  B  I  1 , .25 


HEIGHT=5000.  (FT) 


MACH=.B 


BODY  AXIS 

DATA  -LATERAL- 

N  L 

Y 

BETA 

B. 23792738 

-39.92845659 

- .37444266 

P 

-.00831482 

-2.1 1380615 

.00029122 

R 

-.70740908 

.07970246 

. 00380601 

RUDDER 

-5.0546382S 

9.79342839 

.05309116 

AILERONS 

-1.89279772 

-55.74144715 

.00463283 

CANNARD 

4.38821795 

2.60330899 

.02501018 

THE  A-MATRIX 

IS  : 

-.374442BB 

3S.S2845B59 

3.23792738 

0 . 00000000 

.03182987 
-2. 11380S15 
-.00631482 
1.00000000 

-.99948753 

.07970246 

-.70740908 

.03184556 

.04885258 
0 . 00000000 
0.00000000 
0.00000000 

THE  STATE  VECTOR  TRANSPOSE 
(  BETA  ,  P  ,  R  ,  PHI  ) 

IS  : 

THE  B-MATRIX  IS  '. 

.053091  IS 
9.79342839 
-5.054S382B 
0.00000000 

.004B32S3 

-55.74144715 

-1.89279772 

0 . 00000000 

.02501018 

2.60330899 

4.38821795 

0.00000000 

THE  CONTROL  VECTOR  TRANSPOSE  IS: 
(RUDDER,  AILERON,  CANNARD) 

NEW  FLIGHT  CONDI  TON? 

YES=1.  N0=0.  :  0 

STOP 

024B00  MAXIMUM  EXECUTION  FL. 

0.103  CP  SECONDS  EXECUTION  TIME. 
COMMAND- 


( VEL . 


74 


100  = 
110= 
120= 
'30= 
140= 


PROGRAM  DIMEN( INPUT, OUTPUT, TAPE5=INPUT,TAPEB=0UTPUT) 

REAL  B1 (8) ,B2(8) ,C1 (9) ,C2(9) ,C3(9) ,D1 (9) ,02(9) ,D3(9) 

REAL  El (9) ,E2(9) ,E3(9) ,F1 (9) ,F2(9) ,F3(9) ,G1 (2) ,G2(2) 

REAL  C4 ( 9 ) , C5 ( 9 ) , CB ( 9 ) , 04 ( 9 ) , 05 ( 9 ) , 08 ( 9 ) , F4 ( 9 ) , F5 ( 9 ) , FS ( 9 ) 
REAL  G4 ( 9 ) , G5 ( 9 ) , GB ( 9 ) , AMTL (4,4), BMTL (4,3) 


150  = 

REAL  G3(2) 

,AMAT(4,4) ,BMAT(4,3) 

1S0  =  C 

170= 

PRINT*,"  " 

180= 

PRINT*, " 

THIS  PROGRAM  CONVERTS  NON-DIMENSIONAL- 

190  = 

PRINT*, " 

DATA  TO  DIMENSIONAL-BODY  AXIS-DATA 

FOR 

200= 

500  PRINT*,"  " 

210  = 

PRINT*, " 

SELECT  MODE.  LONGITUDNAL  ONLY  (1)" 

220  = 

PRINT*, " 

LATERAL  ONLY  (2)" 

230  = 

PRINT*, " 

BOTH  ( 3 )  : " 

240= 

READ*, 16 

250=C 

2G0=C 

THE  FIXED  AIRCRAFT  PARAMETERS  ARE  CONTAINED 

IN 

270  =  C 

VECTOR 

B1 . 

280=C 

290  = 

PRINT*, " 

STANDARD  AFTI  AIRCRAFT?" 

300  = 

PRINT*, " 

YES=i.  No=o. 

310  = 

READ*, 11 

320  = 

CALL  AFTI(B1,I1) 

330=C 

340  = 

PRINT*,"  " 

350  = 

PRINT*, " 

FIXED  AIRCRAFT  PARAMETERS" 

3S0= 

PRINT*,"  " 

370= 

PRINT*, " 

WEIGHT  (LBS)  :",B1(1) 

J80= 

PRINT*," 

SPAN  (FT)  :",B1(2) 

•J90  = 

PRINT*, " 

AREA  (FT  **2) :",B1(3) 

400= 

PRINT*, " 

MAC  (FT)  :",B1(4) 

410  = 

PRINT*, " 

I-X  :",B1(5) 

420  = 

PRINT*, " 

I-Y  :",B1(B) 

430  = 

PRINT*, " 

I-Z  :",B1(7) 

440  = 

PRINT*, " 

I-XZ  :",B1(8) 

450= 
480  =  C 
470=C 
480  =  C 
490  =  0 
500=C 
510  = 
520  = 
530  = 
540  = 
550  = 
580  = 
570  = 
580  = 
590  = 
G00  = 
G10= 
B20= 
S30= 


PRINT*, 


THE  PARAMETERS  THAT  VARY  WITH  FLIGHT  CONOITIONS  ARE 
CONTAINEO  IN  THE  VECTOR  B2 

PRINT*,"  " 

PRINT*,"  INPUT  VARING  PARMETERS" 

PRINT*,"  " 

PRINT*,"  HEIGHT  (FT)  :" 

REA0*,B2< 1 ) 

PRINT*,"  MACH 
REA0*,B2(2) 

PRINT*,"  SONIC  VELOCITY  (FT/SEC)  :" 

REA0*,B2(3) 

PRINT*,"  OENSITY  (SLUG/FT  **3)  I" 

REA0*,B2(4) 

PRINT*,"  ANGLE  OF  ATTACK  (OEG)  I" 

REA0*,B2(5) 
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640= 
650= 
660= 
670  = 
680= 
690= 
700  = 
710= 
720= 
730=C 
740=C 
750=C 
760=C 
770=C 
780=C 
790=C 
800=C 
810=C 
820=C 
830=C 
840=6 
850=6 
860  =  6 
670=6 
880=6 
890= 
900=6 
910= 


940  = 
950  = 
960  = 
970  = 
980  = 
990  = 
1000  = 
1010= 
1020= 
1030= 
1040  = 
1050= 
1060  = 
1070  = 
1080  = 
1090* 
1100= 
1110= 
1120= 
1130= 
1140  = 
1150= 
1160=6 
1170=6 
/  180  = 


^920 

(**330 


PRINT*,"  " 

PRINT*,"  MOMENT  60EF." 

PRINT*,"  " 

PRINT*,"  6L  :" 

READ*,B2(6) 

PRINT*,"  6M  :" 

READ*, 82 (7) 

PRINT*,"  6D 
READ*, 82 (8) 

THE  NON  DIMENSI0NAL-5TA8ILITY  AKI9-DATA  ARE  OONTAINED 
IN  THE  OEOTORS  61,  62,  63,  64,  65,  66. 

LATERAL-DIRE6TI0NAL 

61  60NTAINS  THE  N-MOMENT  60EFI6IENTS 

62  60NTAINS  THE  L-MOMENT  60EFI6IENTS 

63  60NTAINS  THE  Y-F0R6E  60EFI6IENT5 

LONGITUDINAL 

64  60NTAINS  THE  D-F0R6E  60EFI6IENT5 

65  60NTAINS  THE  M-MOMENT  60EFI6EINT5 

66  60NTAIN3  THE  L-F0R6E  60EFI6EINT5 

IF<I6.EQ.l)  GO  TO  1200 
PRINT*,"  " 

PRINT*,"  INPUT  THE  LATERAL  STABILITY  GOEFIGIENTS" 
PRINT*,"  " 

PRINT*,"  GN  ,  GL  ,  GY" 

PRINT*,"  BETA  (DEG) 

READ*, 61 ( 1 ) ,62(1 ) ,63(1 ) 

PRINT*,"  P  (RAD) 

READ*, 61 (2) ,62(2) ,63(2) 

PRINT*,"  R  (RAD)  :" 

READ*, 61 (3) ,62(3) ,63(3) 

PRINT*,"  RUDDER  (DEG)  I" 

READ*,G1 (4) ,62(4) ,63(4) 

PRINT*,"  FLAPERON  (DEG)  I" 

READ*, 615, 625, 635 
PRINT*,"  DIFFAIL  (DEG) 

READ*,G15A,G25A,G35A 
PRINT*,"  GANNARD  (DEG) 

READ*,61 (6) ,62(6) ,63(6) 

PRINT*,"  INPUT  ALERON  GOMBINATION  TERMS" 

PRINT*,"  ALERON=  A*FLAPERON  +  B*DIFF  TAIL" 

PRINT*,"  ENTER  A,  B  :" 

READ*,X7,Xa 

61 (5)=(X7*G15)+(X8*G15A) 

62 ( 5 ) = ( X7*G25 ) + ( Xa*G25A ) 

63 ( 5 ) = ( X7*G35 ) + ( X8*G35A ) 


IF(I6.Ea.2)  GO  TO  1300 
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II  H 


1190=C 

1200=  1200  PRINT*, 

1210=  PRINT*,"  INPUT  THE  LONGITUDINAL  STABILITY  COEFICEINTS" 

,';:-.20=  PRINT*,"  " 

1230=  PRINT*,"  CD  CM  CL" 

1240=  PRINT*,"  " 

1250=  PRINT*,"  OELOCITY  (FT/SEC)  I" 

1260=  READ*,C4< 1 > ,C5( 1) ,CB< 1) 

1270=  PRINT*,"  ALPHA  <DEG>  1" 

1280=  READ*,C4(2) ,C5(2) ,CG(2) 

1290=  PRINT*,"  ELEVATOR  (DEG)  I" 

1300=  READ*,C4(3) ,C5<3) ,CG(3) 

1310=  PRINT*,"  T.E.  FLAP  (DEG) 

1320=  READ* , C4 ( 4 ) , C5 ( 4 ) , CS ( 4 ) 

1330=  PRINT*,"  L.E.  FLAP  (DEG)  :" 

1340=  READ*,C4(5) ,C5(5) ,CS(5) 

1350=  PRINT*,"  Q  (RAD)  **  CM  AND  CL  ONLY  ** 

1360=  READ*,C5(6) ,C6(6) 

1370=  PRINT*,"  ALPHA  DOT  (RAD)  **  CM  AND  CL  ONLY  **  I" 

1390=  READ*,C5(7) ,CG(7) 

1390=C 

1400=C  CONVERT  B2(5)  ANGLE  OF  ATTACK  FROM 

1410=C  DEGREES  TO  RADIANS. 

1420=C 

1430=  1300  PI=3. 141592654 
1440=  DTR=180.0/PI 

1450=C 

1460=  B2(5)=B2(5)/DTR 

^70=C 

^80=  CA  =  C0S(B2(5)  ) 

1490=  SA=SIN(B2(5) ) 

1500=C 

1510=C  CONVERT  MOMENTS  OF  INERTIA  FROM  BODY  AXIS  TO 

1520=C  STABILITY  AXIS  TO  MAKE  FUTURE  CALCULATIONS  CONSISTANT. 

1530=C 

1540=  B15=(B1 (5)*(CA**2) ) - ( 2 . 0*B1 ( 8 ) *CA*SA ) + ( B1 (7)*(SA**2) ) 

1550=  B17=(B1 (7)*(CA**2) ) + ( 2 . 0*B1 ( 8 ) *CA*SA ) + ( B1 ( 5 ) * ( SA**2 ) ) 

1560=  B18=( (B1 (5)-Bl (7) ) *CA*SA ) + ( B1 ( 8 ) * ( ( CA**2 ) - ( SA**2 ) ) ) 

1570=C 

1580=  B1(5)=B15 

1590=  B1(7)=B17 

1600=  B1(8)=B18 

1G10=C 

1620=C  CONVERT  Bl(l)  -WE IGHT ( LBS ) -  TO 

1630=C  -MASS(LBS*(SEC**2)/FT>- 

1S40=C 

1650=  Bl(l)=Bl(l)/32.2 


1660=C 

1670=C 

CALCULATE  VELOCITY  -V-  =B2(2)*B2(3) 

1680=C 

B2(2)  -MACH- 

1690=C 

B2(3)  -SONIC  VELOCITY- 

1700=C 

1710= 

V=B2(2)*B2(3) 

1720=C 

-■.30  = 

UO=V*CA 
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1740=  W0=0*SA 

1750=C 

17B0=C  CONVERT  STABILITY  DERIVATIVES  WITH  (1/DEG)  TO  (1/RAD) 

:770=C  SPECIFICALLY  BETA , RUDDER , AILERONS , AND  CANNARDS  LATERAL 

.780=C  ALPHA, ELEVATOR,  AND  FLAPS  LONGITUDNAL 

1790=C 

1800=C 

1810=  IF(IB.EQ.l)  GO  TO  1210 

1820=C 

1830=  Cl ( 1 )=C1 ( 1 )*DTR 

1840=  C2(1)=C2(1)*DTR 

1850=  C3(1)=C3(1)*DTR 

1880=  Cl (4)=C1(4)*DTR 

1870=  C2(4)=C2(4)*DTR 

1880=  C3(4)=C3(4)*DTR 

1880=  Cl (5)=C1 (5)*DTR 

1800=  C2(5)=C2(5)*DTR 

1810=  C3(5)=C3(5)#DTR 

1820=  C1(B)=C1(B)*DTR 

1830=  C2(B)=C2(B)#DTR 

1340=  C3(B)=C3(B)*DTR 

1950=C 

1980=  IF(IB.Ea.2)  GO  TO  1310 

1870=C 

1980=  1210  C4(2)=C4(2)*DTR 
1980=  C5(2)=C5(2)*DTR 

2000=  C6(2)=CB(2)*DTR 

2010=  C4(3)=C4(3)#DTR 

.2020=  C5(3)=C5(3)*DTR 

(1*030=  C8(3)=CS(3)*DTR 

2040=  C4(4)=C4(4)*DTR 

2050=  C5(4)»C5(4)*DTR 

2080=  C8(4)=C8<4)*DTR 

2070=  C4(5)=C4(5)*DTR 

2080=  C5(5)=C5(5)*DTR 

2090=  CB(5)=CB(5)*DTR 

2100=C 
2110=C 

2120=C  G1  G2  G3  CONTAIN  THE  CONVERSION 

2130=C  FACTORS  TO  CONVERT  NON-DIMEN  DATA 

2140=C  TO  DIMEN  DATA  IN  THE  STABILITY  AXIS 

2150=C 

2180=  1310  IFdB.EQ.l)  GO  TO  1220 
2170=C 

2180=  G1 ( 1 )=(B2(4)*B1 (3)*(V*#2)#B1 (2) )/ (2.#B1 (7) ) 

2190=  G1(2)=(B2(4>*B1(3)*V*(B1(2)**2) )/(4.*Bl(7) ) 

2200=  G2( 1 )=G1 ( 1 )*B1 (7)/Bl (5) 

2210=  G2(2)=G1(2)*B1(7)/B1(5) 

2220=  G3(1)=(B2(4)*B1(3)#V)/(2.»B1(1) ) 

2230=  G3(2)=(B2(4)*B1 (3)*B1 (2) )/(4.*Bl ( 1 ) ) 

2240=  CALL  DIMN(C1 ,G1 ,01 ) 

2250=  CALL  OIMN ( C2 , G2 , 02 ) 

2280=  CALL  DIMN ( C3 , G3 , D3 ) 

2270=C 

..7280=C  the  vectors  D1  ,  D2,  D3  CONTAIN  THE 
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220O=C 
2300=C 
2310  = 
320  =  C 
2330=C 
2340=C 
2350=C 
23B0=C 
2370=C 
2380=C 
2390=C 
2400=C 
2410= 
2420  = 
2430  = 


2440  = 
2450  = 
24B0=C 
2470= 
2480  = 
2490  = 
2500  = 
2510= 
2520=C 
2530= 


2540  = 
2550  = 
25B0  = 
^570  = 
^-580=C 


2590  = 
2300  = 
2810= 
2820  = 
2830  = 
2840  = 
2850  = 
2880  = 
2870= 
2880= 
2890  = 
2700  = 
2710= 
2720  = 
2730  = 
2740  = 
2750  = 
2780  = 
2770  = 
2780=C 
2790=C 
2800»C 
2810=C 


2820  = 
830* 


DIMEN  DATA  IN  THE  9TABILITY  AXI3 

CALL  PRM<D1 ,D2,E1 ,E2.B1 ) 

THE  VECTORS  E1,E2  CONTAIN  THE  DIMEN 
DATA  IN  THE  STABILITY  N  AXIS  EXCEPT 
THE  DATA  IS  PRIMED. (  SEE  MC  RUER ,  AIRCRAFT 
DYNAMICS  AND  AUTOMATIC  CONTROL.  PG  257) 

THE  DIMEN  DATA  -BODY  AXIS-  ARE  CONTAINED  IN 
THE  VECTORS  F1,F2.F3 

F1(1)  =  <E1(1) *CA )  +  <  E2 ( 1 ) *SA ) 

F1<2)=(E1 <2)*(CA**2) )-( (El(3>-E2(2) ) #CA*SA ) - ( E2 < 3 ) * ( SA** 

FI (3)=(E1 (3)*(CA**2) )+( (E2(3)+E1 (2) ) *CA*SA ) + ( E2 ( 2 ) * ( SA** 

DO  10  1=4.8 

10  F1(I)=(E1(I)*CA)+(E2(I)*SA) 

F2( 1 )=(E2( 1 )*CA)-(E1 < 1 )*SA) 

F2(2)=(E2(2)*<CA**2) ) - ( < E2 ( 3 ) +E 1 ( 2 ) ) «CA*SA ) + ( El ( 3 ) * ( SA*#2 ) ) 
F2(3)=<E2(3)*(CA»*2) )-( (El(3)-E2(2) ) *CA«SA ) - ( El ( 2 ) * ( SA**2 ) ) 
DO  20  1=4.8 

20  F2(I)=(E2<I)*CA)-(E1<I)*SA> 

F3(1)=D3(1) 

F3(2)=(D3<2)*CA)-(D3(3)*SA) 

F3 ( 3 ) = ( D3 ( 3 ) *CA ) + ( D3 ( 2 ) *SA ) 

DO  30  1=4,8 
30  F3(I)=D3(I) 

PRINT*."  " 

PRINT*,"  HEIGHTS". B2(l) ."  (FT)  MACH=".B2(2) 

PRINT#,"  " 

PRINT*,"  BODY  AXIS  DATA  -LATERAL-" 

PRINT*,"  " 

PRINT*."  N  L  Y" 

PRINT*."  " 

PRINT  50.F1 ( 1 ) ,F2(1 ) .F3( 1 ) 

PRINT  51 .FI (2) .F2(2) .F3(2) 

PRINT  52.F1 (3) .F2(3) .F3{3) 

PRINT  53.F1 (4) .F2(4) ,F3(4) 

PRINT  54.F1 (5) .F2(5) ,F3(5) 

PRINT  55.F1 (8) .F2(6) .F3(8) 

50  FORMAT (5X,"  BETA  " , 5X , 3 ( FI 8 . 9 . 3X ) . "  (VEL.)",/) 

51  F0RMAT(5X,"  P  " . 5X , 3 ( FIB . 9 . 3X ) . / ) 

52  FORMAT (5X."  R  " , 5X , 3 ( F 18 . 9 . 3X ) . / ) 

53  F0RMAT(5X. "RUDDER" ,5X.3(F1B.9,3X) . / ) 

54  F0RMAT(5X. "ALERONS" .5X,3(F1B.9.3X) ./) 

55  F0RMAT(5X. "CANNARD" .5X.3(F1B.9.3X) ) 

PUT  EQUATIONS  OF  MOTION  IN  MATRIX  FORM. 

SEE  MCRUER  PAGE  259 

DO  80  1=1,4 
DO  SO  J  =  1 .4 
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2840=  BO  AMATd,  J)=0.0 
2350=  DO  70  1=1,4 

2860=  DO  70  J=l,3 

2£.-^=  70  BMAT(  I  ,  J)=0.0 

28B0=C 

2880=  U0=0*CA 

2900=  WO=V*SA 

2910=C 

2920=  AMAT(l,l)=F3(l) 

2930=  AMAT(l,2)=(F3(2)+W0)/0  , 

2940=  AMAT( l,3)=<F3<3)-U0)/0 

2950=  AMAT(1,4)=32.174*CA/V 

2980=  AMAT(2,1)=F2(1) 

2970=  AMAT(2r2)=F2(2) 

2980=  AMAT(2,3)=F2(3) 

2990=  AMAT(3,1)=F1(1) 

3000=  AMAT(3,2)=F1 (2) 

3010=  AMAT{3,3)=F1 (3) 

3020=  AMAT(4,2)=1 .0 

3030=  AMAT(4,3)=SA/CA 

3040=  BMAT(1,1)=F3(4) 

3050=  BMAT( 1 ,2)=F3(5) 

3060=  BMAT( 1 ,3)=F3(6) 

3070=  BMAT(2,1)=F2<4) 

3080=  BMAT(2,2)=F2(5) 

3090=  BMAT(2,3)=F2(S) 

3100=  BMAT(3, 1 )=F1 (4) 

3110=  BMAT(3,2)=F1(5) 

3/^=  BMAT(3,3)=F1<G) 

3>oO*C 

3140=  PRINT*,"  " 

3150=  PRINT*,"  THE  A-MATRIX  19  :" 

3180=  CALL  PRINTM<4,4,AMAT) 

3170=  PRINT*,"  " 

3180=  PRINT*,"  THE  STATE  VECTOR  TRANSPOSE  IS  I 

3190=  PRINT*,"  (  BETA  ,  P  ,  R  ,  PHI  )" 

3200=  PRINT*,"  " 

3210=  PRINT*,"  THE  B-MATRIX  IS  I" 

3220=  CALL  PRINTtK  4 , 3 , BMAT ) 

3230=  PRINT*,"  " 

3240=  PRINT*,"  THE  CONTROL  VECTOR  TRANSPOSE  IS 

3250=  PRINT*,"  (RUDDER,  ALERON,  CANNARD)" 

3260=C 

3270=  IF(IS.EQ.2)  GOTO  1320 

3280*C 

3290=C  THE  VECTORS  G4,G5,GS  CONTAIN  THE  FACTORS 

3300=C  THAT  CONVERT  NON  DIMENSIONAL  DATA  TO 

3310=C  DIMENSIONAL  DATA  -STABILITY  AXIS- 

3320=C 

3330=  1220  G4( 1 ) =B2(4)*B1 (3)*V/B1 ( 1 ) 

3340=  G4<2)=G4( 1 )/2. 

3350=  G4 ( 3 ) =B2  <  4 ) *B 1 ( 3 ) * ( V**2 ) / ( 2 . *B 1 ( 1 ) ) 

3360>C 

3570=  Q5 ( 1 ) =B2 ( 4 ) *B 1 ( 3 ) *V*B 1 ( 4 ) /B 1 ( 8 ) 

:-:-:-;-o»  g5<2)=g5<  i  )/2. 


3390  = 
3400= 
3410= 

.;':420=C 
3430= 
3440  = 
3450= 
3460  = 
3470= 
3480=C 
3490=C 
3500=C 
3510=C 
3520=C 
3530=C 
3540=C 
3550=C 
35B0=C 
3570=C 
3580=C 
3590=C 
3600=C 
3610=C 
3620= 
3630= 
3640= 
3650  = 
3660=C 

1^70  = 

Vo80= 
3690=C 
3700  = 
3710= 
3720  = 
3730= 
3740=C 
3750= 
3760  = 
3770=C 
3780= 
3790= 
3800= 
3810= 
3820=C 
3830  = 
3840  = 
3850=C 
3860=C 
3870=C 
3880=C 
3890= 
3900= 
3910= 
3920 »C 
::33o» 


G5(3)=B2(4)*B1 <3)*(B1 (4)**2 ) / ( 4 . #B1 ( 6) > 

G5<4>=G5(3>*V 

G5(5)=B2(4)*B1 (3)*(V**2)*B1 (4)/(2.*Bl (6) ) 

G6(1)=G4(1> 

G6<2)=G4(2) 

GB(3)=B2(4)*B1 <3)*B1 (4)/(4.*Bl ( 1 ) ) 

G6<4)=G6(3)«0 

G6(5)=B2(4)*B1 (3)*<0**2)/(2.*Bl ( 1 ) ) 

THE  VECTORS  D4,D5,DS  CONTAIN  THE  DIMENSIONAL 
DATA  -STABILITY  AXIS- 

*#  NOTE  THE  VECTOR  SEQUENCE  CHANGE  ** 

1  —  VELOCITY  — 

2  —  W  ~ 

3  —  W  DOT  “ 

4  ~  Q  — 

5  —  ELEVATOR  — 

6  —  T.E.  FLAP  - 

7  —  L.E.  FLAP  — 

D4(1)=G4(1)*<-B2(8)-C4(1)  ) 

D4<2)=G4(2)*(B2(B)-C4(2) ) 

D4(3)=0.0 

D4(4)=0.0 

DO  700  I=5r7 

700  D4(I)=G4<3)*<-C4(I-2) ) 

D5( 1 )=G5< 1 )*(B2(7)+C5< 1 ) ) 

D5(2)=G5(2)*C5(2) 

D5(3)=G5<3)*C5<7) 

D5(4)=G5(4)*C5<S) 

DO  710  I=5»7 
710  D5( I )=G5(5)*C5( 1-2) 

DS(1)=GG(1)*<-B2<6)-CG(1)  ) 

DS ( 2 ) =G6 ( 2 ) *  < -C6 ( 2 ) -B2 ( 8 ) ) 

D6(3)=G6(3)*CS(7) 

D6<4)=G6<4)«(-CG(6) ) 

DO  720  I=5r7 

720  D6(I)=G6(5)*(-CB(I-2)) 

THE  VECTORS  F4,F5rF6  CONTAIN  THE 
DATA  IN  BODY  AXIS 

F4( 1 )=(D4( 1 )*(CA**2) )-( (D4(2)+DB{ 1 ) ) *SA*CA ) + ( D6 ( 2 ) * ( SA** 
F4(2)=(D4(2)*{CA**2) )+( (D4( 1 )-DG(2) ) *SA*CA ) - ( DS < 1 )*(SA** 
F4 ( 3 ) * ( D4 ( 3 ) *  <  CA**2 ) ) - ( DS ( 3  >  *SA*CA ) 

DO  800  1=4.7 
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3940=  800  F4( I )=(D4< I )*CA)-(DB( I )*SA) 

3950=C 

3960=  F5 ( 1 )  =  <  D5 ( 1 ) *CA ) - <  D5 ( 2 ) *SA ) 

•  370=  F5 ( 2 )  =  ( D5 ( 2 ) *CA )  +  ( D5 ( 1 >  *SA ) 

j980=  F5(3)=D5(3)*CA 

3990=C 

4000=  DO  810  1=4,7 

4010=  810  F5(I)=D5(I) 

4020=C 

4030=  FB(1 )=(DS< 1 )*(CA**2) )-< (D6(2)-D4( 1) ) *SA*CA ) - ( D4( 2 ) * ( SA**2 ) ) 

4040=  F6<2)={D6(2)*(CA**2) )+( (D6( 1 )+D4(2> ) *SA*CA ) + ( D4 ( 1 )*(SA**2) ) 

4050=  FB(3)=(D6(3)*(CA**2) )+(D4(3)*SA*CA) 

4060 =C 

4070=  DO  820  1=4,7 

4080=  820  FB(I)=<D6(I)*CA)+(D4(I)*SA) 

4090  =  C 

4100=  PRINT*,"  " 

4110=  PRINT*,"  HEIGHT  I", 82(1),"  MACH  I", 82(2) 

4120=  PRINT*,"  " 

4130=  PRINT*,"  80DY  AXI8  DATA  -LONGITUDINAL-" 

4140=  PRINT*,"  " 

4150=  PRINT*,"  X  M  Z" 

4160=  PRINT*,"  " 

4170=  PRINT  90,F4( 1 ) ,F5( 1 ) ,FB( 1 ) 

4180=  PRINT  91 ,F4(2) ,F5(2) ,F6(2) 

4190=  PRINT  92, F4(3) ,F5(3) ,FB(3) 

4200=  PRINT  93, F4(4) ,F5(4) ,F6(4) 

4210=  PRINT  94, F4(5) ,F5(5) ,F6<5) 

^20=  PRINT  95,F4(6)  ,F5(B)  ,F6<6) 

1^30=  PRINT  96, F4(7)  ,F5(7)  ,F6(7) 

4240=C 

4250=  90  F0RMAT(5X, "VELOCITY" ,3X,3(F16.9,3X) ,/) 

4260=  91  F0RMAT<5X,"  U  " , 3X , 3(F1B . 9 , 3X ) , / ) 

4270=  92  FORMAT (5X,"  W-DOT  " , 3X , 3 ( FIB . 9 , 3X ) , / ) 

4280=  93  FORMAT (5X,"  Q  " , 3X , 3 ( FIB . 9 , 3X ) , / ) 

4290=  94  F0RMAT(5X, "ELEVATOR", 3X,3<F1B. 9, 3X) ,/) 

4300=  95  F0RMAT(5X, "T.E.FLAP" ,3X,3(F1B.9,3X) ,/) 

4310=  96  F0RMAT<5X,"L.E.FLAP",3X,3(F1B.9,3X) ,/) 

4320=C 

4330=  23=1.0/ ( 1.0-FG(3) ) 

4340=C 

4350=  AMTL (1,1) =F4 ( 1 ) + ( F4 ( 3 ) *FB ( 1 ) *F6 ( 3 ) ) 

4360=  AMTL  < 1 , 2 )  =  <  F4  <  2 ) *U0 )  +  ( F4 ( 3 ) *U0*F6 ( 2 ) *Z3 ) 

4370=  AMTL ( 1 , 3 )  =  <  F4  <  4 ) -WO )  +  ( F4 ( 3 ) * ( F6 ( 4 ) +U0 ) *Z3 ) 

4380=  AMTL(1,4)=-32.174*CA*(1.0-(F4(3)*Z3) ) 

4390=C 

4400=  AMTL(2,1)=F6(1)*Z3/U0 

4410=  AMTL(2,2)=F6(2)*Z3 

4420=  AMTL(2,3)=(F6(4)+U0)*23/U0 

4430=  AMTL(2,4)=-32. 174*CA*Z3/U0 

4440=C 

4450*  AMTL (3,1) =F5 ( 1 ) + ( F5 ( 3 ) *F6 ( 1 ) *Z3 ) 

4460=  AMTL ( 3 , 2 ) = ( F5 ( 2 ) *U0 ) + ( F5 ( ^ ) *U0*F6 ( 2 ) *Z3 ) 

4  470=  AMTL ( 3 , 3 ) =F5 ( 4 )  +  ( F5 ( 3 ) * ( F6 ( 4 ) +U0 ) *Z3 ) 

:;-j480»  AMTL(3,4)=-32. 174*CA*F5(3)*Z3 

4490*0 
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4500  = 

4510  = 

4520  = 

’  530  = 

■ 54o=C 
4550= 

45G0  = 
4570= 
4580=C 
4590  = 
4600= 

4610  = 

4620=C 

4630= 

4640  = 

4650  = 
46G0=C 
4670  = 

4680  = 

4680  = 
4700=C 
4710=C 
4720=  1320 
4730  = 

4740  = 

4750  = 

4760  = 
4770=C 
4780 =C 
^V90=C 


AMTL<4, 1 ) =0.0 
AMTL(4,2)=0.0 
AMTL(4,3)=1 .0 
AMTL<4,4)=0.0 

BMTL( 1 , 1 )=F4(5)+<F4(3)*FG(5)*Z3) 
BMTLd  ,2)=F4(6)  +  <F4(3)*FG<6)*Z3) 
BMTL( 1 ,3)=F4(7)+(F4<3)*F6(7)#Z3) 

BMTL(2, 1 )=FG(5)*Z3/U0 
BMTL(2,2)=F6(6)*Z3/U0 
BMTL ( 2  r  3 ) =FB ( 7 ) *Z3/U0 

BMTL(3, 1 ) =F5(5)+(F5(3)*F6(5)*Z3) 
BMTL ( 3 , 2 ) =F5 ( 6 ) + ( F5 ( 3 ) *FG ( 6 ) *Z3 ) 
BMTL(3,3)=F5(7)+<F5(3)#F6<7)»Z3) 

BMTL (4, 1 ) =0.0 
BMTL(4,2)=0.0 
BMTL(4,3)=0.0 


PRINT*,"  " 

PRINT*,"  NEW  FLIGHT 
PRINT*,"  YES=1.  NO 
READ*, 12 

rF<I2.EQ.l)  GO  TO  500 


NEW  FLIGHT  CQNDITON?" 
YEs=i.  NO=o. 


4800  = 

STOP 

4810  = 

END 

4820=C 
4830  = 

SUBROUTINE  D IMN ( A 1 , G 1 , A2 ) 

4840  = 

REAL  A1 (9) ,G1 (2) ,A2(g) 

4850= 

A2(1)=A1(1)*G1(1) 

4860= 

DO  10  J=2,3 

4870  = 

10 

A2(J)=A1( J)*G1(2) 

4880  = 

DO  20  K=4,G 

4890  = 

20 

A2(K)=A1 <K )*G1 ( 1 ) 

4800  = 

RETURN 

4810  = 

END 

4920=C 
4930  = 

SUBROUTINE  PRM  ( A1  ,  A2  , El  ,  E: 

4840  = 

REAL  A1(9),A2<9),E1(9),E2 

4850  = 

Pl=l-( (B1 (8)**2) / (B1 <7)*B 

4980  = 

P2=B1 (8)/Bl (5) 

4970  = 

P3=B1 (8)/BJ (7) 

4980  = 

DO  10  1=1,6 

4990  = 

El ( I )  =  ( A1 ( I )  +  <  P3*A2< I ) ) )/l 

5000  = 

10 

E2( I ) = (A2( I )+( P2*A1 ( I ) ) )/| 

5010  = 

RETURN 

5020  = 

END 

5030=C 
•040  = 

SUBROUTINE  AFTI(B1,I1) 
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5050  = 

REAL  81(7) 

5oeo= 

IF(Il.EQ.l)  GO  TO  100 

5070=C 

;:)80= 

PRINT*,"  " 

5030  = 

PRINT*,"  INPUT'FIXED  AIRCRAFT  PARAMETERS 

5100= 

PRINT*,"  " 

5110  = 

PRINT*,"  WEIGHT  (LBS)  " 

5120  = 

READ*, 81 ( 1 ) 

5130  = 

PRINT*,"  " 

5140  = 

PRINT*,"  SPAN  (FT)  I" 

5150  = 

READ*,B1 (2) 

5160= 

PRINT*,"  " 

5170  = 

PRINT*,"  AREA  (FT  **2) 

5130  = 

READ*,B1 (3) 

5130  = 

PRINT*,"  " 

5200= 

PRINT*,"  MEAN  AERODYNAMIC  CHORD  (FT)  I" 

5210  = 

READ*,B1 (4) 

5220  = 

PRINT*,"  " 

5230= 

PRINT*,"  MOMENTS  OF  INERTIA  " 

5240  = 

PRINT*,"  " 

5250  = 

PRINT*,"  I-X  :" 

5260  = 

READ*,B1(5) 

5270  = 

PRINT*,"  " 

5280  = 

PRINT*,"  I-Y 

5230= 

READ*,B1 (6) 

5300  = 

PRINT*,"  " 

5310  = 

PRINT*,"  i-z 

5320  = 

READ*,B1 (7) 

^30= 

PRINT*,"  " 

*340  = 

PRINT*,"  i-xz 

5350  = 

READ*,B1 (8) 

5360  = 

GO  TO  200 

5370  = 

100 

81  ( 1 )=21018.0 

5380  = 

B1 (2) =30.0 

5330  = 

B1 (3) =300.0 

5400= 

B1 (4)=11 .32 

5410  = 

Bl(5)=10033.423 

5420  = 

B1 (6)=5387S.2SS 

5430  = 

B1 (7)=61278.452 

5440= 

Bl(8)=282. 13217 

5450  = 

200 

RETURN 

5460  = 

END 

5470=C 

5480  = 

SUBROUTINE  PRINTM ( N , L , A ) 

5430  = 

REAL  A(4,4) 

5500= 

PRINT*,"  " 

5510=C 

5520= 

DO  100  I=1,N 

5530  = 

PRINT  200,(A(I,J),J=1,L) 

5540  = 

100 

CONTINUE 

5550=C 

5560= 

200 

FORMAT( 1X,8(F15.8,3X) ) 

5570  = 

RETURN 

5580  = 

END 

530=*E0R 
:5600  =  *E0F 
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Appendix  B-1 

Copg>uter  Program  to  calculate  System  Zeros 
As  part  of  this  thesis  effort.  Program  ZERO  is  developed  to  calculate 
the  zeros  of  the  linear  matrix  system 


i:--. 


X  =  A  X  +  B  11 
2.  =  C  X  +  D  U 


(B.l) 


J  cr 


'.'•J 


*1 


ZERO  calculates  transmission  auid  decoupling  zeros  by  solving  the  appropriate 
generalized  eigenvalue  problem  with  the  powerful  QZ  technique  (Ref  6) . 

Three  specific  eigenvalue  problems  are  used  (Ref  8). 

The  generalized  eigenvalue  problem  to  calculate  trauismission  zeros 
depends  on  the  number  of  inputs,  m,  and  the  number  of  outputs,  l,  of  the 
system. 

For  f  ^  m  the  problem  is 


(B.2) 


'a  B  0  ■ 

t- 

o 

o 

H 

n 

£  D  0 

z  =  A 

0  0  0 

O 

O 

1 _ 

o 

o 

o 

_ 1 

For  m  ^  f  the  problem  is 

A  B 

S.  E 


m 


0 

0 


z  =  X 


0 

0 


o 

o 


(B.3) 


O 


0  0 


Die  generalized  eigenvalue  problem  to  calculate  decoupling  zeros  is 


"l  0  O 

n 

z  «  X  0  0  0 

.0  0  0 

In  the  above  generalized  eigenvalue  problems,  k  is  a  randcxti  m  -  by  -  i 
real  matrix. 

'Hie  generalized  eigenvalue  problem  can  be  described  as 

z  -  X  z  (B.5) 

The  QZ  algorithm  has  three  steps.  The  first  reduces  the  matrix  M  to 
upper  Hessenberg  form  wd  the  matrix  ^  to  upper  triwgular  form.  Next 
the  Hessenberg  matrix  is  further  reduced  to  quasi- triangular  form,  while 
maintaining  the  triangulau:  form  of  the  other  matrix.  The  final  step 
reduces  the  quasi-triemgular  matrix  further,  such  that  the  remaining 
2  -  by  -  2  blocks  correspond  to  pairs  of  conplex  eigenvalues,  nie  QZ 
algorithm  returns  qucuitities  whose  ratio  gives  the  generalized  eigenvalues. 

The  conputer  program  ZERO  makes  scxne  decisions  before  the  system 
zeros  are  finally  claissified.  Refer  to  the  flow  chart  shown  in  Figure 
B-1.  For  each  path  of  the  zero  progreun,  the  appropriate  generalized 
eigenvalue  problem,  G.E.P.,  is  solved  twice  (runs  1  emd  2)  with  only  the 
K  matrix  changing.  The  outputs  of  runs  1  and  2  are  conqpeured.  If  the 
value  of  em  eigenvalue  is  the  s^une  for  both  runs,  then  it  is  classified 
as  a  zero. 
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The  zeros  calculated  by  the  QZ  algorithm,  set  z, 


may  contain  some 


TA' 

decoi^ling  zeros.  Hence,  set  z  is  compared  with  the  decoupling  zeros 

X  A 

calculated  by  the  QZ  algorithm,  set  z^^.  Those  values  in  set  z^^  that 

also  exist  in  set  z^  are  eliminated  from  set  z_. .  The  rem£d.ning  zeros  in 

D  TA 

set  z^^  form  the  set  of  treuismission  zeros,  set  z^. 

Esqperience  shows  that  occasionally  some  values  2u;e  incorrectly 
identified  as  system  zeros.  Hence,  ZERO  uses  three  different  tests  in 
order  to  filter  out  incorrect  values  from  the  sets  z^  and  z^,  and  to 
further  classify  the  real  decot^ling  zeros. 

Test  A  (Subroutine  Syszero)  uses  singular  value  decooposition  to 

0 

check  the  rank  of  the  system  matrix 


P(X)  = 


XI  -  A,  -B 
n  —  — 


(B.6) 


An  invariemt  zero  causes  P(X)  to  lose  rank.  All  treuismission  zeros  eure 
invariwt  zeros,  along  with  some  decoupling  zeros.  All  of  the  reeil 
tremsmission  zeros  and  real  decoupling  zeros  are  sent  through  Test  A. 
Any  treuismission  zero  that  does  not  cause  P^(X)  to  lose  reuik  is  rejected 
6is  a  zero.  Any  decoupling  zero  that  does  caiise  P^(X)  to  lose  rank  is 
classified  as  eui  invariant  decoupling  zero. 

Test  B  (Subroutine  BDZ)  uses  singular  value  deccm^osition  to  check 


the  rank  of  the  matrices 


All  of  the  real 


decoi^ling  zeros  are  sent  through  Test  B,  in  order  to  classify  them  as 
input,  output,  or  input/output  decoupling  zeros.  Input  decot^ling  zeros 
cause  the  matrix  fx I  -A,b1  to  be  rauik  deficient.  Output  decoupling  zeros 


cause  the  matrix  I  n  —  I  to  be  ramk  deficient.  Input/output  decoupling 


zeros  cause  both  matrices  to  be  rank  deficient.  A  real  decoupling  zero 
that  catises  neither  matrix  to  lose  remk  is  rejected  as  a  zero. 

Test  C  (Subroutine  Complex)  is  for  complex  system  zeros.  All  conqplex 
zeros  must  occur  in  conjugate  ped.rs.  Any  conplex  zero  without  a  conjugate 
is  rejected. 

ZERO  is  aui  interactive  program.  The  required  inputs  are  listed  in 

Tedile  B-1.  Input  data  may  be  read  from  Tape!  =  ZDATA.  The  order  of  the 

inputs  for  Tape3  ZDATA  is  eJLso  listed  in  Tad)le  B-1.  The  ZDATA  input 
2 

format  is  (312,  (n  -l-njiH-2ni)F20.7) .  Two  ASD  libraries,  EISPACK  emd  IHSL, 
must  be  connected  before  executicxi  of  the  program  ZERO. 

The  current  version  of  ZERO  is  not  100%  accurate.  Sometimes  an 
incorrect  value  will  pass  through  the  filter  tests.  These  values  can 
be  eeisily  identified  because  they  are  very  large  values.  In  the  third 
exaii^>le  the  number  -19,722,484.728  is  incorrectly  identified  as  a 
trsuismisslon  zero. 


TABLE  B-1 


Required  Input  Data  for  Program  Zero 


Symbol 


Dimension 


Input  Order 
Interactive  ZDATA 


Number  of  States 
A  Matrix 
Number  of  Inputs 
B  Matrix 

Number  of  Outputs 
C  Matrix 
D  Matrix 
F  Matrix 


n-by-n 


n-by-m 


£-by-n 

£-by-m 

£,-by-n 


THIS  TEST  PROGRAM  CALCULATES  ZEROS  OF  A  LINEAR 
MATRIX  SYSTEM. 

MAXIMUM  #  OF  STATES  =  8. 

OF  INPUTS  =  4. 

OF  OUTPUTS=  4. 

TO  INCREASE  MAXIMUM,  SEE  INSTRUCTIONS  AT  THE 
BEGINNING  OF  THE  PROGRAM  LISTING. 

FILE  ZEROLIS, ID=LEWIS,SN=AFIT 

ARE  THE  MATRICES  STORED  ON  TAPE  3  =  ZDATA  ? 

(  1  =  YES.  0  =  NO  >  :  0 

ENTER  THE  DIMENSION  OF  THE  SYSTEM  :  6 

ENTER  THE  A  MATRIX  <B  BY  G)  I 

ROW  1  :  1 , 0 , 0 , 0 , 0 , 0 
ROW  2  :  0,1, Or 0,0,0 
ROW  3  :  0,0, 3, 0,0,0 
ROW  4  :  0,0, 0,-4, 0,0 
ROW  5  :  0,0,0, 0,-1 ,0 
ROW  6  :  0,0, 0,0, 0,3 

ENTER  THE  NUMBER  OF  INPUTS  I  2 

ENTER  THE  B  MATRIX  (S  BY  2) 

ROW  1  :  0,-1 
ROW  2  :  -1,0 
ROW  3  :  1,-1 
ROW  4  :  0,0 
ROW  5  :  0,1 
ROW  G  :  -1,-1 

ENTER  THE  NUMBER  OF  OUTPUTS  I  3 

ENTER  THE  C  MATRIX  (3  BY  G) 

ROW  1  :  1,0,0, 1,0,0 
ROW  2  :  0,1, 0,1, 0,1 
ROW  3  :  0,0, 1,0,0, 1 

IS  THERE  A  D  MATRIX  ? 

(YES=1  N0=0)  :  0 

THE  D  MATRIX  IS  A  (3  BY  2)  ZERO  MATRIX 
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HE  A  MATRIX  IS 


1.0000 

0 . 0000 

0.0000 

0.0000 

0.0000 

0 . 0000 

0.0000 

1.0000 

0.0000 

0.0000 

0 . 0000 

0 . 0000 

0.0000 

0.0000 

3.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-4.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-1.0000 

0.0000 

0.0000 

0.0000 

0 . 0000 

0.0000 

0.0000 

3.0000 

THE  B 

MATRIX  IS  : 

0.0000 

- 1 . 0000 

-1.0000 

0 . 0000 

1 . 0000 

-1.0000 

0.0000 

0.0000 

0.0000 

1.0000 

-1.0000 

-1.0000 

THE  C  MATRIX  IS  : 

1.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0 . 0000 

0.0000 

1.0000 

0.0000 

1.0000 

0.0000 

1 . 0000 

0.0000 

0.0000 

1 .0000 

0.0000 

0 . 0000 

1 . 0000 

***  WARNING 

ERROR 

(lER  = 

33)  FROM 

IMSL  ROUTINE 

LSVDF 

**■*  WARNING 

ERROR 

(lER  = 

33)  FROM 

IMSL  ROUTINE 

LSVDF 

WARNING  33  MEANS  A  MATRIX  DOES  NOT 

HAVE  FULL  RANK.  THIS  MESSAGE  IS  EXPECTED 


THE  TRANSMISSION  ZEROS 

REAL  PART  IMMAGINARY  PART 

2.0000000000  0. 

THE  DECOUPLING  ZEROS 

-1.0000  OUTPUT  INVARIANT 

-4.0000  INPUT 

DO  YOU  WISH  TO  CHANGE  MATRICES 
AND  CONTINUE? 
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STOP  =  0 

CHANGE  A  ONLY  =  ’l 
CHANGE  B  ONLY  =  2 
CHANGE  C  ONLY  =  3 
CHANGE  ArBrC,0  =4 
CHANGE  B,C  ONLY  =  5 
CHANGE  C,D  ONLY  =  6 
4 

ENTER  THE  DIMENSION 

ENTER  THE  A  MATRIX 

ROW  1  :  0,1, 0,0, 0,0 
ROW  2  :  0,0, 1,0, 0,0 
ROW  3  :  0,0, 0,0, 0,0 
ROW  4  :  0,0, 0,0, 1,0 
ROW  5  :  0,0, 0,0, 0,1 
ROW  6  :  0,0, 0,0, 0,0 

ENTER  THE  NUMBER  OF 

ENTER  THE  B  MATRIX 

ROW  1  :  0,0 
ROW  2  :  0,0 
ROW  3  :  1,0 
ROW  4  :  0,0 
ROW  5  :  0,0 
ROW  6  :  0,1 

ENTER  THE  NUMBER  OF 

ENTER  THE  C  MATRIX 

ROW  1  :  1,1, 0,0, 0,0 
ROW  2  :  0,0,0,-l ,-l 

’IS  THERE  A  D  MATRIX 
(YES=1  N0=0)  :  1 

ENTER  THE  D  MATRIX 

ROW  1  :  1,0 
ROW  2  :  1,0 


OF  THE  SYSTEM  :  6 
(6  BY  6)  : 


INPUTS  :  2 
(6  BY  2) 

OUTPUTS  :  2 
2  BY  S) 

0 

? 

2  BY  2) 
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THE  A 

MATRIX  IS  : 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0 . 0000 

0.0000 

1.0000 

0.0000 

0.0000 

0 . 0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0 . 0000 

0.0000 

0.0000 

0 . 0000 

0.0000 

1.0000 

0.0000 

0 . 0000 

0.0000 

0 . 0000 

0.0000 

0 . 0000 

1.0000 

0 . 0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

THE  B 

MATRIX  IS  : 

0 . 0000 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

THE  C 

MATRIX  IS  : 

1.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-1.0000 

0.0000 

-1.0000 

0.0000 

0.0000 

THE  D 

MATRIX  IS  : 

1.0000 
1 . 0000 

0.0000 

0.0000 

WARNING  ERROR 
***  WARNING  ERROR 

WARNING  33  MEANS 
HAVE  FULL  RANK. 

(lER  = 
<IER  = 
A  MATRIX  DOES 
THIS  MESSAGE 

33)  FROM 
33)  FROM 
NOT 

IS  EXPECTED 

IMSL  ROUTINE 
IMSL  ROUTINE 

LSVDF 

LSVDF 

THE  TRANSMISSION  ZEROS 

REAL  PART  IMMAGINARY  PART 


.34116390191 
.3411S390191 
l.OOOOOOC  00 
.  682327803 


1 . 1815414000 
-1.1615414000 

0. 

0. 


THERE  ARE  NQ  DECOUPLING  ZEROS 
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THIS  TEST  PROGRAM  CALCULATES  ZEROS  OF  A  LINEAR 
MATRIX  SYSTEM. 

MAXIMUM  #  OF  STATES  =  8. 

OF  INPUTS  =  4. 

OF  OUTPUTS=  4. 

TO  INCREASE  MAXIMUM,  SEE  INSTRUCTIONS  AT  THE 
BEGINNING  OF  THE  PROGRAM  LISTING. 

FILE  ZEROLIS, ID=LEWIS,SN=AFIT 

ARE  THE  MATRICES  STORED  ON  TAPE  3  =  ZDATA  ? 

(  1  =  YES.  0  =  NO  )  :  1 

THE  A  MATRIX  IS  : 


0.0000  0.0000  1.0000  .0318 

.0489  -.3744  .0318  -.9995 

0.0000  -39.9283  -2.1138  .0797 

0.0000  8.2379  -.0083  -.7074 


THE  B  MATRIX  IS  ! 


0.0000  0.0000  0.0000 

.0531  .0048  .0250 

9.7929  -55.7388  2.8033 

-5.0544  -1.8927  4.3880 


THE  C  MATRIX  IS  : 


0.0000  1.0000  0.0000  0.0000 

1.0000  0.0000  0.0000  0.0000 

0.0000  0.0000  0.0000  1.0000 

***  WARNING  ERROR  (lER  =  33)  FROM  IMSL  ROUTINE  LSODF 

***  WARNING  ERROR  (lER  =  33)  FROM  IMSL  ROUTINE  LSUDF 

WARNING  33  MEANS  A  MATRIX  DOES  NOT 
HAVE  FULL  RANK.  THIS  MESSAGE  IS  EXPECTED 


THE  TRANSMISSION  ZEROS 

REAL  PART  IMMAGINARY  PART 


-19722484.728  0. 

THERE  ARE  NO  DECOUPLING  ZEROS 


I 
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w 


0  YOU  WISH  TO  CHANGE  MATRICES 
AND  CONTINUE? 

STOP  =  0 

CHANGE  A  ONLY  =  1 
CHANGE  B  ONLY  =  2 
CHANGE  C  ONLY  =  3 
CHANGE  A,B,C,D  =  4 
CHANGE  B,C  ONLY  =  5 
CHANGE  C,D  ONLY  =  B 
:  3 

CHANGE  C-MATRIX  TO  F-MATRIK  ? 
YES=1,  N0=0  :  1 


THE  C  MATRIX  IS  ! 


0.0000 
1.0000 
0 . 0000 


1 . 0000 
0.0000 
0.0000 


0 . 0000 
.2500 
0.0000 


0.0000 
.OOSO 
1 . 0000 


•»#*  WARNING  ERROR  (lER  =  33)  FROM  IMSL  ROUTINE  LSODF 

WARNING  33  MEANS  A  MATRIX  DOES  NOT 
HAVE  FULL  RANK.  THIS  MESSAGE  IS  EXPECTED 


THE  TRANSMISSION  ZEROS 


REAL  PART 


IMMAGINARY  PART 


-4.0000000000 


THERE  ARE  NO  DECOUPLING  ZEROS 


DO  YOU  WISH  TO  CHANGE  MATRICES 
AND  CONTINUE? 


STOP  =  0 
CHANGE  A  ONLY 
CHANGE  B  ONLY 
CHANGE  C  ONLY 
CHANGE  A,B,C,D 
CHANGE  B,C  ONLY 
CHANGE  CrD  ONLY 
0 
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100= 

110= 

120= 

..30= 

140= 

150=C 

160= 

170=C 

180=C 

190=C 

200=C 

210=C 

220=C 

230=C 

240=C 

250=C 


260=C 
270=C 
280=C 
2S0=C 
300=C 
310=C 
320=C 
330=C 
340=C 
350=C 
360=C 
370  =  C 


^ao*c 

^90=C 

400=C 


410  = 


420  = 


430  = 

440  = 

450  = 

460= 

470  = 

480= 

490= 

500= 

510= 

520= 

530= 

540= 

550=C 

560= 

570= 


PROGRAM  2ER0( INPUT, OUTPUT, TAPE5=INPUT,TAPE6=aUTPUTr 
1  ZDATA,TAPE3=ZDATA) 

REAL  TRAR( 16) ,TRAI (16) ,OECR< 16) ,OECI ( 16) ,FMAT<8,8) 

REAL  AMAT (8,8), BMAT (8,8), CMAT (8,8), DMAT (8,8) 

INTEGER  IFLTdG)  ,IFL(16), 10(16) 

LOGICAL  MATO 

THI8  PROGRAM  CAN  BE  SIMPLY  ALTERED  TO  ACCOMMIDATE 
MORE  STATES  AND  INPUTS/OUTPUTS. 

THIS  IS  DONE  BY  CHANGING  DIMENSION  STATEMENTS 
THE  IMPORTANT  DIMENSION  NUMBERS  ARE! 

N  -  NUMBER  OF  STATES 
M  -  NUMBER  OF  INPUTS 
L  -  NUMBER  OF  OUTPUTS 
NM  -  N  +  MAX(M,L) 

NN  -  N  +  M  +  L 

CHANGE  ALL  ARRAYS  I4ITH  (8,8)  TO  (N,N) 

CHANGE  ALL  ARRAYS  WITH  (12,12)  TO  (NM,NM) 

CHANGE  ,12,  IN  IMSL  FUNCTION  LSODF  IN  SUBROTINE  RDZ 
AND  SUBROUTINE  SYSZERO  TO  ,NM. 

CHANGE  ALL  VECTORS  WITH  (18)  AND  ARRAYS 
WITH  (16,16)  TO  (NN)  OR  (NN,NN) 

ALSO  CHANGE  THE  ,16,  IN  THE  QZ  SUBROT I NES 
TO  ,NN. 

A  SAMPLE  COMMAND  IN  EDITOR  IS  I  /16/=/NN/ , A,U,0  (VETO) 
SAVE  NEW  FILE  THEN  RECOMPLILE 


PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
READ*,N1 
IF(N1 .EG. 1 ) 


THIS  TEST  PROGRAM  CALCULATES  ZEROS  OF  A  LINEAR" 
MATRIX  SYSTEM." 

MAXIMUM  #  OF  STATES  =  8." 

OF  INPUTS  =4." 

OF  OUTPUTS=  4." 

TO  INCREASE  MAXIMUM,  SEE  INSTRUCTIONS  AT  THE" 
BEGINNING  OF  THE  PROGRAM  LISTING." 

FILE  ZEROLIS, ID=LEWIS,SN=AFIT" 


ARE  THE  MATRICES  STORED  ON  TAPE  3 
(  1  =  YES.  0  =  NO  )  :  " 

GO  TO  1100 


ZDATA  ?" 


PRINT*,"  " 
MATCH=0 


580=  5400  PRINT*,"  ENTER  THE  DIMENSION  OF  THE  SYSTEM  :" 
590=  READ*,N 


GOO=C 

S10=C  INPUT  THE  BLOCK  A  MATRIX 

620=C 


-<530=  PRINT*,"  " 
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640*  PRINT*, “  ENTER  THE  A  MATRIX  (",N,"  BY  ",N,")  I” 

650*  CALL  READM(NrN,AMAT) 

160*  IF(MATCH.EQ. 1 )G0  TO  5430 

B70*C 

680*  PRINT#,"  " 

690*  5410  PRINT*,"  ENTER  THE  NUMBER  OF  INPUTS  I" 

700*  READ*,M 

710*C 

720»C  INPUT  THE  BLOCK  B  MATRIX 

730=C 

740*  PRINT*,"  " 

750*  PRINT*,"  ENTER  THE  B  MATRIX  {",N,"  BY  ",M,")" 

760*  CALL  READM<N,M,BMAT) 

770*  IF ( MATCH. EQ. 2) GO  TO  5430 

780*C 

790*  PRINT*,"  " 

800*  5420  IF(Nl.EQ.l)  GO  TO  5421 

810*  5424  PRINT*,"  ENTER  THE  NUMBER  OF  OUTPUTS  I" 

820=  READ*,L 

830*  GO  TO  5422 

840*  5421  PRINT*,"  CHANGE  C-MATRIX  TO  F-MATRIX  ?" 

850*  PRINT*,"  YES=1,  N0=0  I" 

860*  READ*,N2 

870*  IF<N2.EQ.l)  GO  TO  5423 

880*  GO  TO  5424 

890*  5422  CONTINUE 

900*C 

|,>aiO*C  INPUT  THE  BLOCK  C  MATRIX 

li!20*C 

930*  PRINT*,"  " 

940=  PRINT*,"  ENTER  THE  C  MATRIX  (",L,"  BY  ",N,")" 

950*  CALL  READM(L,N,CMAT) 

960*C 

970=  IF (MATCH. EG. 4)  GO  TO  5435 

980*  IF(MATCH.Ea.6)  GO  TO  5435 

990*  IF ( MATCH. NE.O)  GO  TO  1200 

1000*C 

1010*  5435  PRINT*,"  " 

1020*  PRINT*,"  IS  THERE  A  D  MATRIX  ?" 

1030*  PRINT*,"  (YES*1  N0*0)  I" 

1040*  READ*, 19 

1050*  IFdS.EQ.l)  GO  TO  5440 

1060*C 

1070*  DO  5450  I=1,L 

1080*  DO  5450  J*1,M 

1090*  5450  DMATd,  J)*0.0 
1100*  PRINT*,"  " 

1110*  PRINT*,"  THE  D  MATRIX  IS  A  (",L,"  BY  ",M,")  ZERO  MATRIX" 

1120*  PRINT*,"  " 

1130*  GO  TO  5430 

1140*C 

1150*  5440  PRINT*,"  " 

1160*  PRINT*,"  ENTER  THE  D  MATRIX  (",L,"  BY  ",M,")" 

■I;  70*  CALL  READM(L,M,DMAT) 
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1180=C 

490=  5430  GO  TO  1200 
.200=0 
1210=0 

1220=0  READ  THE  INPUT  DATA  FROM  TAPE  3 

1230=0 

1240=0 

1250=  1100  OONTINUE 
1260=0 

1270=  READ(3,45)  N»M,L 

1280=  45  FORMAT (312) 

1280=0 

1300=  DO  50  I=1,N 

1310=  DO  50  J=lrN 

1320=  50  READ{3,55)  AMAT(IrJ) 

1330=0 

1340=  DO  60  1=1, N 

1350=  DO  60  J=1,M 

1360=  60  READ(3,55)  BMAT(I,J) 

1370=0 

1380=  DO  70  I=1,L 

1380=  DO  70  J=1,N 

1400=  70  READ(3,55)  OMAT(I,J) 

1410=0 

1420=  DO  80  1=1, L 

1430=  DO  80  J=1,N 

M40=  80  READ<3,55)  FMAT<I,J) 

[¥50=0 

1460=  19=0 

1470=  DO  80  1=1, L 

1480=  DO  90  J=1,M 

1490=  90  DMATd,  J)=0.0 

1500=0 

1510=0 

1520=  55  FORMAT (F20. 7) 

1530=  GO  TO  1200 

1540=0 

1550=  5423  OONTINUE 

1560=  DO  65  1=1, L 

1570=  DO  85  J=1,N 

1580=  85  OMAT(I,J)=FMAT(I,J) 

1590=0 

1600=  1200  OONTINUE 
1610=0 
1620=0 

1630=  IF(MAT0H.Ea.2)  GO  TO  91 

1640=  IF(MATGH.Ea.3)  GO  TO  92 

1650=  IF(MAT0H.EQ.5)  GO  TO  91 

1660=  IF(MATGH.EQ.B)  GO  TO  92 

1670=0 

1680=  PRINT*,"  " 

1690=  PRINT*,”  THE  A  MATRIX  IS  I" 

.  :00=  PRINT*,"  " 

1/10=  OALL  PRINTM(N,N,AMAT) 


ZDATA 
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1720=C 

1730=  IF(MATCH.EQ. 1 )  GO  TO  95 

:  ,;40=C 

V;50=  91  PRINT*,"  " 

1760=  PRINT*,"  " 

1770=  PRINT*,"  THE  B  MATRIX  19  I" 

1780=  PRINT*,"  " 

1790=  CALL  PRINTM(N,M,BMAT) 

1800=C 

1810=  IF(MATCH.EQ.2)  GO  TO  95 

1820=C 

1830=  92  PRINT*,"  " 

1840=  PRINT*,"  " 

1850=  PRINT*,"  THE  C  MATRIX  IS  I" 

1860=  PRINT*,"  " 

1870=  CALL  PRINTM(L,N,CMAT) 

1880=  PRINT*,"  " 

1880=  IF (MATCH. EG. 3)  GO  TO  95 

1900=  IF(MATCH.EQ.5)  GO  TO  95 

1910=C 

1920=  IF(19.EQ.l)  GO  TO  96 

1930=  GO  TO  95 

1940=  96  PRINT*,"  " 

1950=  PRINT*,"  THE  D  MATRIX  IS  I" 

1960=  PRINT*,"  " 

1970=  CALL  PRINTM(L,M,DMAT) 

1980=  PRINT*,"  " 

^0=  95  CONTINUE 

l-foo=c 

2010=C 

2020=C  CALCULATE  ZEROS 

2030=C 

2040=C 

2050=  NN=N+M+L 

2060=  MATV=. FALSE. 

2070=C 

2080=C 

2090*C  SUBROUTINE  TZ  CALCULATES  THE  TRANSMISSION  ZEROS 

2100=C 

2110=C 

2120=  CALL  TZ<AMAT,BMAT,CMAT,DMAT,NN,N,M,L,MATO, 

2130=  +TRAR,TRAI , IFLT) 

2140*C 

2150=C  DOUBLE  CHECK  REAL  TRANSMISSION  ZEROS  BY  TESTING 

21S0=C  THE  RANK  OF  THE  SYSTEM  MATRIX.  A  TRANSMISSION 

2170=C  ZERO  CAUES  THE  SYSTEM  MATRIX  TO  BE  RANK  DEFICIENT. 

2180=C 

2190=  EX3*1.0E-3 

2200»C 

2210*  DO  440  1=1, NN 

2220=  IF( IFLT( I ) .EQ.O)GO  TO  440 

2230=  IF<TRAI( I ) .LT.EX3.AND.TRAI(I> .GT.-EX3)  GO  TO  445 

2240=  GO  TO  440 

*50*  445  Z»TRAR(I) 
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22B0*  CALL  SYSZERO( AMAT,BMAT,CMAT,DMAT,Z»N,M,L, IT, 1 ) 

2270*  IF(IT.EQ.O)  IFLT(I)»0 

:2B0*  440  CONTINUE 

2290*C 

2300*C  SUBROUTINE  OZ  CALCULATES  THE  DECOUPLING  ZEROS 

2310*C 

2320*C 

2330*  CALL  DZ(AMAT.BMAT.CMAT.DMAT.NNf N,M.LrNATU, 

2340*  -t-DECR^OECIrlFL) 

2350*C 

23S0*C 

2370*C  COMPARE  THE  TRANSMISSION  ZEROS  FOUND  IN  TZ 

2380*C  WITH  THE  DECOUPLING  ZEROS  FOUND  IN  DZ. 

23S0=C  IF  THE  SAME  VALUE  APPEARS  IT  IS  ONLY  A 

2400*C  DECOUPLING  ZERO. 

2410*C 

2420*  EX*1.0E-5 

2430=C 

2440*  DO  200  I*1,NN 

2450*C 

24B0*  IF(  IFLT(  I )  .EQ.O)  GO  TO  200 

Z470*  ATRAR*ABS<TRAR< I) ) 

2480=  ATRAI*ABS(TRAI(I) ) 

2490=C 

2500*  DO  210  J=1,NN 

2510=C 

2520=  IF(IFL( J) .EQ.O)  GO  TO  210 

1^30=  AOECR= ABS  ( DECR  (  J  )  ) 

^rei40=  AOECI*ABS(DECI(J) ) 

2550=C 

25B0»C  COMPARE  REAL  SYSTEM  ZEROS 

2570=C 

2580*  IF(ATRAI .GT.EX3)  GO  TO  201 

2590=  IF<ADECI.GT.EX3)  GO  TO  210 

2S00=  IF(ATRAR.LT.EX3)  GO  TO  202 

2B10  =  C 

2B20=  DT* ( TRAR ( I ) -DECR ( J ) ) /TRAR ( I ) 

2B30*  IF(ABS(DT) .LT.EX)  IFLT(I)*0 

2B40*  GO  TO  210 

2850*C 

2BB0*  202  IF(TRAR(I).LT.DECR<J)+EX.AND.TRAR(I).GT.DECR(J)-EX)  IFLT(I)=0 

2B70*  GO  TO  210 

2680*C 

2G90=C  COMPARE  IMMANGINARY  SYSTEM  ZEROS 

2700*C 

2710*  201  IF(ADECI.LT.EX3)  GO  TO  210 

2720*  IF(ATRAR.LT.EX3)  GO  TO  203 

2730*C 

2740*  DT I  * ( TRAR ( I ) -DECR ( J ) ) /TRAR ( I ) 

2750*  IF(ABS(DTI) .LT.EX)  GO  TO  204 

2760-  GO  TO  210 

2770-C 

2780*  203  IF(TRAR(I)  .LT.DECR(  J)-*-EX.AND.TRAR<I)  .GT.DECRl  J)-EX)  GO  TO  204 

,90*  GO  TO  210 
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Z800=C 

>310= 

204 

2620= 

2830= 

2840= 

2850=C 

2860  = 

205 

2870=C 

2880= 

210 

2890= 

200 

2S00»C 

2910=C 

2920=0 

2930=C 

2940=C 

2950= 

2960  = 

2970=C 

2980= 

2990  = 

3000= 
3010  = 
3020= 
3030= 
3040= 
3050= 


60= 

70=C 


3080= 

3090= 


3100= 

3110= 

3120= 

3130  = 

3140= 

3150  = 

3160  = 

3170=C 

3180=C 

3igo=c 

3200»C 

3210=C 

3220*C 

3230=C 

3240= 

3250= 


3260= 
3270= 
3280  = 
3290= 
3300= 
3310  = 

:-2o» 

3330= 


305 

300 


315 

310 


550 


IF(ATRAI.LT.EX3)  GO  TO  205 
DTII=(TRAI(I)-DECI(J) >/TRAI(I) 

IF<ABS(DTII) .LT.EX)  IFLT(I)=0 
GO  TO  210 

IF ( TRA I ( I ) . LT . DEC I ( J ) +EX . AND . TR A I ( I ) . GT . DEC I ( J ) -EX ) 

CONTINUE 

CONTINUE 

CHECK  CALCULATED  COMPLEX  ZEROS  FOR  COMPLEX 
CONJIGATE  PAIRS.  COMPLEX  ZEROS  MUST  OCCUR  IN 
CONJIGATE  PAIRS. 

CALL  COMPLEX(NN,TRAI,IFLT) 

CALL  COMPLEX (NN> DEC I r IFL) 


KTZ=0 

ITZ=1 

DO  300  1=1 fNN 

IF< IFLT( I ) .EG. 1 )  GO  TO  300 
KTZ=KTZ+1 

IF(KTZ.EQ.NN)  GO  TO  305 

GO  TO  300 

ITZ=0 

CONTINUE 

KDZ=0 
IDZ  =  1 

DO  310  J=1,NN 

IF( IFL( J) .EG. 1 )G0  TO  310 

KDZ=KDZ+1 

IF(KDZ.EO.NN)  GO  TO  315 

GO  TO  310 

IDZ=0 

CONTINUE 


TEST  DECOUPLING  ZEROS  TO  DETERM IN  WETHER  THEY 
ARE  INVARIANT,  AND  IF  THEY  ARE  AN  INPUT  OR  OUTPUT 
ZERO. 


DO  500  I=1,NN 
10<l)=0 

IF( IFL( I ) .EG.O)  GO  TO  500 

IF(DECI ( I ) .LT.EX3.AND.DECI ( I ) .GT.-EX3)  GO  TO  550 

GO  TO  500 

Z=DECR(I) 

CALL  RDZ ( AMAT , BMAT , CMAT ,Z,N,M,L,Kia) 

CALL  SYSZERO{AMAT,BMAT,CMAT,DMAT,Z,N,M,L, ID,2) 
IF(ID.EG.l)  GO  TO  510 
IF<KI0.E0.1)  I0(I)=1 


IFLT( I )»0 
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3340= 

3350= 

r^360= 

:  >70= 

3380= 

3390= 

3400= 

3410=C 

3420= 

3430= 

3440= 

3450  = 

3480  = 

3470  = 

3480=C 

3490=C 

3500= 

3510= 

3520  = 

3530= 

3540= 

3550  = 

35S0=C 

3570  = 

3580  = 

3590= 

3600  = 

3610= 

^0=C 

OT30= 

3640  = 

3650  = 

3660=C 

3670= 

3680= 

3690= 

3700= 

3710= 

3720  = 

3730= 

3740= 

3750=C 

3760= 

3770= 

3780* 

3790  = 

3800= 

3810  = 

3820  = 

3830= 

3840= 

3850= 

3880= 

3R70« 


IF<KI0.EG.2)  I0(I)=2 
IF(KI0.EQ.3)  I0<I)=3 
GO  TO  500 

510  IF<KI0.EQ.1>  I0(n=4 
IF(KI0.EQ.2)  I0(I)=5 
IF(KI0.EQ.3>  I0<I)=6 
500  CONTINUE 

IF( ITZ.EQ.O.AND. IDZ.EQ.O)  GO  TO  399 


WARNING  33  MEAN8  A  MATRIX  DOES  NOT" 

HAVE  FULL  RANK.  THIS  MESSAGE  IS  EXPECTED" 


PRINT*, " 
PRINT*, " 
PRINT*," 
PRINT*, " 
399  CONTINUE 


IF(ITZ.EQ.O)  GO  TO  400 
PRINT*,"  " 

PRINT*,"  THE  TRANSMISSION  ZEROS" 
PRINT*,"  " 

PRINT  412 
PRINT*,"  " 

DO  410  1=1, NN 

IFdFLTd)  .EQ.O)  GO  TO  410 

PRINT*,"  " 

PRINT  411,(TRARd),TRAId)) 

410  CONTINUE 

415  IFdDZ.EQ.O)  GO  TO  405 
PRINT*,"  " 

PRINT*,"  THE  DECOUPLING  ZEROS" 


DO  420  I=1,NN 
IFdFLd  )  .EQ.O)  GO  TO  420 
IF( I0< I ) .EQ.O)  GO  TO  421 
IFdOd  )  .EQ.  1  )  GO  TO  422 
IFd0d).Ea.2)  GO  TO  423 
IFdOd). EQ. 4)  GO  TO  424 
IFdOd)  .EQ. 5)  GO  TO  426 
IFdOd). EQ. 6)  GO  TO  427 

PRINT*,"  " 

PRINT  413,DECRd) 

GO  TO  420 
427  PRINT*,"  " 

PRINT  417,DECRd) 

GO  TO  420 
426  PRINT*,"  " 

PRINT  418,DECRd) 

GO  TO  420 
424  PRINT*,"  " 

PRINT  419,DECRd) 

GO  TO  420 
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II  II 


3880=  423  PRINT*, 

3880=  PRINT  414rDECR(I) 

N  :;300=  GO  TO  420 

3810=  422  PRINT*,"  “ 

3320=  PRINT  418,DECR(I) 

3830=  GO  TO  420 

3340=  421  PRINT*,"  " 

3850=  PRINT  412 

3360=  PRINT  411,(DECR(I),DECI(I)) 

3370=  420  CONTINUE 

3380=C 

3330=  411  FORMAT (2G20. 11 ) 

4000=  412  FORMAT (8X, "REAL  PART" , 8X, “ IMMAGINARY  PART") 

4010=  413  F0RMAT(5X,F10.4,5X, "  INPUT/OUTPUT") 

4020=  414  F0RMAT<5X,F10.4,5X, "  OUTPUT") 

4030=  416  F0RMAT(5X,F10.4,5X, "  INPUT") 

4040=  417  F0RMAT(5X,F10.4,5X, "  INPUT/OUTPUT  INUARIANT") 

4050=  418  F0RMAT(5X,F10.4,5X, "  OUTPUT  INVARIANT") 

4080=  413  F0RMAT(5X,F10.4,5X, "  INPUT  INVARIANT") 

4070=C 

4080=C 

4030=C 

4100=  GO  TO  430 

4110=  400  PRINT*,"  " 

4120=  PRINT*,"  THERE  ARE  NO  TRANSMISSION  ZEROS" 

4130=  PRINT*,"  " 

4140=  PRINT*,"  " 

g*m50=  PRINT*,"  " 

VbO=  go  to  415 


4170  = 

405 

PRINT*," 

II 

4180= 

PRINT*, " 

THERE  ARE  NO  DECOUPLING  ZEROS" 

4130= 

PRINT*, " 

II 

4200= 

PRINT*, " 

II 

4210= 

PRINT*, " 

II 

4220= 

4230=C 

4240=C 

430 

CONTINUE 

4250=C 

SELECT 

OPTION.  CONTINUE  OR  RETURN 

42S0=C 

4270= 

PRINT*, " 

II 

4280= 

PRINT*, " 

DO  YOU  MISH  TO  CHANGE  MATRICES 

4230= 

PRINT*, " 

AND  CONTINUE?" 

4300= 

PRINT*, " 

II 

4310= 

PRINT*," 

STOP  =0" 

4320= 

PRINT*, " 

CHANGE  A  ONLY  =  1" 

4330= 

PRINT*," 

CHANGE  B  ONLY  =2" 

4340= 

PRINT*," 

CHANGE  C  ONLY  =  3" 

4350= 

PRINT*, " 

CHANGE  A,B,C,D  =  4" 

4360= 

PRINT*, " 

CHANGE  B,C  ONLY  =5" 

4370= 

PRINT*, " 

CHANGE  C,D  ONLY  =6" 

4380= 

PRINT*," 

•  II 

■ 

4330= 

READ*, MATCH 

4400= 

IF (MATCH 

.EQ.l)  GO  TO  5400 

10= 

IF<MATCH 

.EQ.2)  GO  TO  5410 

$ 


4420* 

4430= 

^  >0* 

^50* 

4460= 

4470= 

4480=C 

4490=C 

4500= 

4510= 

4520= 

4530=C 

4540= 

4550= 

4560= 

4570  = 

4580=C 

4590= 

4600= 

4610=C 

4620= 

4630= 

4640= 

4650=C 

4660= 

4670= 

4680= 


4710= 

4720  = 

4730=C 

4740=C 

4750= 

4760  = 

4770=C 

4780= 

4790= 

4800= 

4810=C 

4820= 

4830  = 

4840=C 

4850=C 

4860= 

4870= 

4880= 

4890=C 

4900*C 

4910=C 

4920-C 

4930>C 

4940= 

.  50= 


IF(MATCH.Ea.3)  GO  TO  5420 
IF(MATCH.EQ.4)  GO  TO  5400 
IF{MATCH.Ea.5)  GO  TO  5410 
IF (MATCH. EG. 6)  GO  TO  5420 
STOP 
END 


SUBROUTINE  READM(L>N.A) 
REAL  A (8, 8) 

PRINT*,"  " 

DO  100  I=lrL 
PRINT*,"  ROW  ",I,“ 
READ*,<A(I,J>,J=1,N) 

100  CONTINUE 

RETURN 

END 

SUBROUTINE  PRINTM ( L , N , A ) 
REAL  A (8, 8) 

PRINT*,"  " 

DO  100  1=1, L 
PRINT  200,(A<I,J),J=1,N) 
100  CONTINUE 

200  F0RMAT(1X,8(F12.4,3X)) 
RETURN 
END 


SUBROUTINE  MSUB < A , B , C , L ,N ) 
REAL  A(8,8),B(8,8),C(8,8) 

DO  100  1  =  1, L 
DO  100  J=1,N 
100  C(I,J)=A(I,J)-B(I, J) 

RETURN 

END 


SUBROUT 1 NE  ABCD ( AA , BB , A , B , C , D , NN , N , M , L  > 

REAL  AA(16,16),BB(16,16),A<8,8),B(8,8), 
-»-C(8,8),D(8,8) 

THIS  SUBROUTINE  SETS  UP  THE  MATRICES  AA  AND  BB 
WHERE  AA=(  A,B  ),  AND  BB=  (  1,0  ) 

(  C,D  ),  (  0,0  ) 

DO  200  I=1,NN 
DO  200  J=1,NN 
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4960=  AA<I,J)=0. 

4970=  200  BB(IrJ)=0. 

-  :\}o=c 

4£rgo=  DO  210  1  =  1, N 

5000=  210  BB<I,I)=1. 

5010=C 

5020=  DO  220  I=1,N 

5030=  DO  220  J=1,N 

5040=  220  AA(I,J)=A<I, J) 

5050=C 

5060=  DO  230  I=1,N 

5070=  DO  230  J=1,M 

5080=  230  AA( I rN+J)=B( I , J) 

5090=C 

5100=  DO  240  1  =  1, L 

5110=  DO  240  J=1,N 

5120=  240  AA(N+I , J)=C( I , J) 

5130=C 

5140=  DO  250  1=1, L 

5150=  DO  250  J=1,M 

5160=  250  AA(N+I,N+J)=D(I,J) 

5170=C 

5190=  RETURN 

5190=  END 

5200=0 
5210=0 

5220=  3UBR0UTINE  aZALGM(N,A,B,ALFR,ALFI ,BETA, 

raO=  +MATZ,Z,IFAIL) 

6Z40=  logioal  matz 

5250=0 

5260=0  THIS  8UBR0UTINE  0ALL3  THE  QZ  SUBROUTINES  FROM  THE 

5270=0  SLATEO  (OR  EISPAOK)  LIBRARY 

5280=0 

5290=  REAL  A < 16 , 16 ) , B ( 16 , 16 ) , ALFR( 16) , ALFI ( 16) , BETA ( 16 ) , Z ( 16 , 16 ) 

5300=  ISAVE=IFAIL 

5310=  IFAIL=1 

5320=  ESP=-1.0E-6 

5330=  OALL  GZHES( 16,N,A,B,MATZ,Z) 

5340=  OALL  GZIT < 16 , N , A , B , ESP , MATZ, Z , IFAIL ) 

5350=  OALL  QZVAL( 16,N, A,B, ALFR,ALFI ,BETA,MATZ,Z) 

5360*  IF( IFAIL. NE.O)RETURN 

5370=  IF(MATZ)OALL  aZVEO( 16,N, A,B, ALFR, ALFI ,BETA,Z) 

5380=  RETURN 

5390=  END 

5400=0 

5410=  SUBROUTINE  OAZERO ( AA , BB , NN ,MM, MATO , ZERR , ZERI , 

5420=  +TR,TI ,LABR,LABI , IFL, JF, IX) 

5430=0 

5440=0 

5450=  REAL  AA( 16, 16) ,BB( 16, 16) ,ALFR( 16) , ALFI ( 16) ,BETA( 16) , 

5460=  +00(16,16) ,LABR( 16) ,LABI ( 16) ,ZERR( 16) ,ZERI(16) , 

5470=  +TR( 16) ,TI( 16) 

5480=  INTEGER  IFL (16) 

r'.-  30=  LOGIOAL  MATO 


5500=C 

5510=C 

r  ’>0=C 

S«JO=C 

5540=C 

5550=C 

55B0=C 

5570=C 

5580=C 

5530=C 

5600=C 

5610-C 

5620=C 

5630=C 

5640=C 

5650=C 

5660- 

5670= 

5660=C 

5690  = 

5700=C 

5710  = 

5720= 

5730  = 

5740  = 

5750= 

5760= 


5790= 

5800= 

5810= 

5820= 

5830= 

5840  = 

5850  = 

5860  = 

5870  = 

5880= 

5890  = 

5900  = 

5910= 

5920= 

5930= 

5940= 

5950= 

5960= 

5970=C 

5980*8 

5990  =  C 

6000  = 

6010= 

6020= 

r  -lO- 


THI3  SUBROUTINE  TAKES  DATA  COMPUTED  BY  THE  QZ 
ALGORITHM  AND  CALCULATES  THE  SYSTEM  ZEROS.  IT 
ALSO  DECIDES  IF  THE  NUMBER  IS  OR  ISN'T  A  SYSTEM 
ZERO.  THIS  IS  DONE  BY  USING  THE  QZ  ALGORITHM 
TWICE  AND  COMPARING  THE  OUTPUTS  OF  EACH  RUN.  IF 
THE  VALUE  IS  SAME  FOR  BOTH  RUNS,  IT  IS  A  SYSTEM 
ZERO. 


THE  VALUE  OF  IX  SIGNALS  WHICH  SYSTEM  ZERO 
IS  BEING  CALCULATED. 

TRANSMISSION  ZEROS  IX=1 
DECOUPLING  ZEROS  IX=2 

IF(IX.EQ.l)  EXl=1.0E-a 
IF<IX.EQ.2)  EXl=1.0E-4 

EX2=1 .OE-8 

IFAIL=1 

CALL  QZALGM ( NN , AA , BB , ALFR , ALF I , BETA , 
+MATV,VV, IFAIL) 

IF( IFAIL.EQ.O)  GOTO  150 
PRINT*,"  RUN  FAILED" 

150  CONTINUE 

DO  120  1=1, NN 
LABR( I )=0. 

IF(JF.EQ.l)  TR<I)=0. 

LABI < I )=0. 

IF(JF.EQ.l)  TI(I)=0. 

BT=BETA(I) 

IFL( I  )=0 

IF(ABS(BT) .GE.EX2)  GOTO  40 
GO  TO  120 
40  IFL(I)=1 

LABR(I)=ALFR<I)/BT 
LABI ( I )=ALFI ( I )/BT 
IF(JF.EQ.2.)  GOTO  120 
TR< I )=LABR(I ) 

TI(I)=LABI(I) 

120  CONTINUE 

IF(JF.EQ.2)  GOTO  50 
GOTO  1000 
50  CONTINUE 

COMPARE  CALCULATED  ZEROS  FROM  RUNS  1  AND  2 

DO  130  I=1,NN 

IF( IFL( I ) .EQ.O)  GOTO  130 

2ERR( I)=LABR( I ) 

ZERI(I)=LABI<I) 


S040=C 

6050  = 

30= 

6070 =C 

6080= 

60S0=C 

6100  = 

6110=C 

6120= 

6130  = 

6140=C 

6150=C 

6160=C 

6170=C 

6180  = 

6130= 

6200= 

6210= 

6220= 

6230=C 

6240= 

310 

6250= 

6260= 

6270=C 

6280 =C 

G230=C 

6300= 

400 

6330  = 

6340=C 

6350  = 

410 

6360  = 

6370= 

6380=C 

6380= 

420 

6400= 

6410= 

6420-C 

6430= 

430 

6440= 

6450=C 

6460= 

43 

6470  = 

42 

6480= 

6480  = 

130 

6500=C 

6510=  : 

1000 

6520  = 

6530  = 

6540=C 

8550=C 

6560= 

:;7o»c 


AZR»ABS(ZERR(I)) 

AZI==ABS<ZERI<I)} 


K=0 


DO  42  J-1,NN 

ATR=ABS(TR( J) ) 
ATI»ABS(TI( J)  ) 


COMPARE  REAL  ZEROS 

IF(AZI.GT.EX2)  GO  TO  400 
IF(ATI.GT.EX2)  GO  TO  43 
IF(A2R.GT.EX2)  GO  TO  310 

IF(ZERR<I).LT.TR(J)+EX1.AND.ZERR(I).GT.TR(J)-EX1>  GO  TO  42 
GO  TO  43 

DT=  <  ZERR ( I ) -TR  <  J ) ) /ZERR ( I > 

IF(ABS<DT) .LT.EXl )  GO  TO  42 
GO  TO  43 

COMPARE  IMMANGINARY  ZEROS 

IF(ATI.LT.EX2)  GO  TO  43 
IF(AZR.GT.EX2)  GO  TO  410 

IF(ZERR(I).LT.TR(J)+EX1.AND.ZERR(I).GT.TR(J)-EX1)  GO  TO  420 
GO  TO  43 

DT=<ZERR(I)-TR( J) )/ZERR(I) 

IF (ABS(DT) .LT.EXl )  GO  TO  420 
GO  TO  43 

IF(AZI .GT.EX2)  GO  TO  430 

IF(2ERI(I).LT.TI(J)+EX1.AND.ZERI(I).GT.TI<J)-EX1)  GO  TO  42 
GO  TO  43 

DTI=(ZERI(I)-TI(J))/ZERI(I) 

IF<ABS(DTI ) .LT.EXl )  GO  TO  42 


K=K  +  1 
CONTINUE 

IF(K.Ea.NN)  IFL(I)=0 
CONTINUE 

CONTINUE 

RETURN 

END 


SUBROUTINE  TZ( A,B,C,DrNNrNrMrLrMATV,TRARrTRAI , IFLT) 


108 


6580*C 

6590»  REAL  AA( 16 , 16 ) ^ BB ( 1 6 > 16 ) , A < 8 r 8 ) r B ( 8 r 8 ) r C ( 8 r 8 ) , D ( 8 > 8 > r 

':  “00=  +TRAR(  16)  ,TRAI  (16)  ,TR(  16)  »TI(  16)  ,LABR(16)  rLABI  (  16 )  ,  W  ( 16 , 1 6  ) 

^j10=  integer  IFLT(i6) 

66Z0=C 

6630-  LOGICAL  MATO 

BG40=C 

6650=C  THIS  SUBROUTINE  AND  THE  SUBROUTINES  CALLED  BY  IT 

SSBO=C  COMPUTE  THE  TRANSMISSION  ZEROS  OF  A  SYSTEM 

S670=C 

6680=  DO  100  11=1,2 

6690=  CALL  ABCD ( AA , BB , A , B , C , D ,NN,N , M , L ) 

6700=  IF(L.LT.M)  GO  TO  200 

S710=C 

6720=  DO  110  I=1,L 

6730=  110  AA(N+IrN+M+I)=l. 

6740=C 

6750=  GO  TO  300 

S760=C 

6770=  200  CONTINUE 

6780=C 

6790=  DO  210  1=1, M 

6800=  210  AA(N+L+I ,N+I )=1 . 

6810=C 

6820=  300  CONTINUE 

6830=C 

6840=  DO  220  I=1,M 

P950=  DO  220  J=1,L 

fjlBO=  220  AA<N+L+1,N+M+J)=RANF(DUM)*10. 

BB70=C 

6880=C 

6890=  JF=II 

6900=  CALL  CAZER0(AA,BB,NN,00,MATV,TRAR,TRAI , 

6910=  +TR,TI ,LABR,LABI , IFLT, JF,1 ) 

6920=  100  CONTINUE 

B930=C 

6940=  RETURN 

6950=  END 

6960=C 

6970=  SUBROUTINE  DZ < A , B , C , D , NN , N , M , L , MATV , DECR , 

6980=  +DECI,IFL) 

6990=C 

7000=C 

7010=  INTEGER  IFLdS) 

7020=C 

7030=C  THIS  SUBROUTINE  AND  THE  SUBROUTINES  CALLED  BY 

7040=C  IT  COMPUTE  THE  DECOUPLING  ZEROS  OF  A  SYSTEM. 

7050=C 

7060=  REAL  AA ( 16 , 16 ) , BB ( 16 , 1 6 ) , A ( 8 , 8 ) , B < 8 , 8 ) , C ( 8 , 8 ) 

7070=  +,D(8,8) ,TR( 16) ,TI ( 16) ,VV( 16, IS) , 

7080=  ZLABR( 16) ,LABI ( 16 ) , DECR ( 16 ) , DECI (16) 

7090=  LOGICAL  MATO 

7100=C 

"10=  DO  100  H  =  l,2 
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7120= 
7130=C 
-40= 
-.50  = 
71B0=C 
7170= 
7180  = 
7190=C 
7200= 
7210= 
7220= 
7230=C 
7240=C 
7250  = 
7280  = 
7270= 
7280  = 
7290  = 
7300= 
7310=C 
7320= 
7330  = 
7340  = 
7350  = 
73B0=C 
7370=C 
73eo=c 
^0=C 
{5)0=0 
^10  =  C 
7420  =  C 
7430=C 
7440  = 
7450= 
74B0=C 
7470= 
7480= 
740O  =  C 
7500= 
7510= 
7520= 
7530=C 
7540  = 
7550= 
75B0=C 
7570  = 
7580=C 
7590= 
7800  = 
7810  = 
7820  = 
7830=0 
7840= 
-750= 


CALL  ABCD ( AA  r  BB r  A , B , C  t  D r NN , N , M , L) 

DO  no  I  =  nL 
110  AA(N+I ,N+M+I )  =  1  . 

DO  120  I  =  nM 
120  AA(N+L+I ,N+I )  =  1  . 

DO  130  I=lrM 
DO  130  J  =  nL 

130  AA ( N+L+ 1 , N+M+ J ) =R ANF ( DUM  >*10. 


JF=II 

CALL  CAZERO ( AA , BB , NN , 00 , MATO , DECR , DEC I , 
+TR,TI ,LABR,LABI , IFL, JF,2> 

100  CONTINUE 
RETURN 
END 


SUBROUT I NE  RDZ (A,B»C,Z,N,M,LrKIO) 

REAL  A<8r8)>B(8,8}.C(8r8)rS<8)tS0(8> 

REAL  U(12>12),U0(12n2),SI(8r8>,SIA(8,8) 
REAL  DCIN( 12. 12) ,DCOT( 12,12) ,WK< 12) ,WKO( 12) 


THIS  SUBROUTINE  DETERMINES  WEATHER  A  DECOUPLING  ZERO 
IS  AN  INPUT, OUTPUT, OR  INPUT/OUTPUT  DECOUPLING 
ZERO.  THIS  IS  DONE  BY  TAKING  THE  SINGULAR  OALUE 
DECOMPOSITION  OF  THE  MATRICES  (SI-A,B)  AND  (SI-A) 

(  C  ) 


TO  TEST  FOR  FULL  RANK 


MO=N+M 

LO=N+L 


EX1=1 .OE-1 
EX2=1 .OE-5 


DO  190  1=1, N 
DO  ISO  J=1,N 
190  SI(I,J)=0. 

DO  200  I=1,N 
200  SI(I,I)=Z 


CALL  MSUB(SI ,A,SIA,N,N) 

DO  210  1=1, N 
DO  210  J=1,N 
DCIN( I , J)=SIA( I  ,  J) 

210  DCOT(I,J)=SIA(I,J) 

DO  220  I=1,N 
DO  220  J=1,M 


no 


76B0= 
7B70=C 
7B80= 
l  ;  30= 
■^00= 
7710=C 
7720= 
7730=C 
7740= 
7750= 
77B0= 
7770= 
7780= 
7790= 
7800= 
7810  = 
7820= 
7830=C 
7840= 
7850=C 
7880= 
7870= 
7880= 
7890= 
7900= 
7910= 
7920= 
7930= 

7980  = 
7970= 
7980= 
7990= 
8000=C 
8010=C 
8020  = 
8030= 
8040= 
8050=C 
8080=C 
8070=C 
8080=C 
8090=C 
8100=C 
8110=C 
8120=C 
8130=C 
8140  = 
8150=C 
8180= 
8170= 
8180=C 
8190= 


220  DCIN(I>J-i-N)=B(Ir  J) 

DO  230  1=1 rL 
DO  230  J=1,N 
230  DCOT(I+N, J)=C<I, J) 

CALL  LSODF(DCIN, 12,N,M0rUf 12f0rS,MK, lER) 

KI0=0 

KI=0 

DO  300  1=1 fN 

IF(S(I) .GT.EXl)  GO  TO  300 
IF<S(I) .LT.-EXl)  GO  TO  300 

IF(S( I ) .LT.EX2.AND.S( I ) .GT.-EX2)  GO  TO  350 
GO  TO  300 
350  KI=1 
300  CONTINUE 

CALL  LSODF  <  DCOT  1 1 2  r  LO  t  N , UO  r 1 2  r  0 , SO  >  MK  0  r I ERO  > 
K0T=0 

DO  310  1  =  1, N 

IF(SO(I) .GT.EXl)  GO  TO  310 
IF(SO( I ) .LT.-EXl )  GO  TO  310 

IF(S0(I).LT.EX2.AND.S0(I).GT.-EX2)  GO  TO  380 
GO  TO  310 
380  K0T=1 
310  CONTINUE 

IF<KI.EQ.l)  KI0=1 
IF<KOT.EQ.l)  KI0=2 
IF(KI.Ea.l.AND.KOT.EQ.l)  KI0=3 
RETURN 
END 


SUBROUTINE  SYSZERO < A , 8 , C ,D , Z,N , M , L , IZ , IX ) 

REAL  A(8,8),B<8,8;rC<8,8)rD(8,8),S(8),MK(12) 

REAL  SYSMdZ,  12)  ,SI  (8,8)  ,SIA(8,8)  ,U(  12, 12) 

THIS  SUBROUTINE  DETERMINS  THE  INOARIANT  ZEROS 
ALL  TRANSMISSION  ZEROS  ARE  INVARIANT 
SOME  DECOUPLING  ZEROS  ARE  INVARIANT 

THE  SUBROUTINE  TAKES  THE  SYSTEM  MATRIX  (  SI-A  ,  -B) 

(  C  ,  D) 

AND  TESTS  FOR  FULL  RANK  BY  SINGULAR  VALUE  DECOMPOSITION 
AN  INVARIANT  ZERO  CAUSES  THE  SYSTEM  MATIX  TO  LOOSE  RANK 

IZ=0 

MO”N+M 

LO-N+L 

IF(IX.EQ.Z)  GO  TO  20 


t 

u 
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8200=C 

8210= 

EX1=1 .OE-1 

8220  = 

EX2=1 .OE-9 

-30= 

GO  TO  25 

6240=C 

8250= 

20 

EX1=1 .OE-1 

8280= 

EX2=1 .OE-5 

8270  = 

25 

CONTINUE 

8280=C 

8280= 

DO  100  I=1,L0 

8300= 

DO  100  J=1,M0 

8310= 

8320=C 

100 

SY3M(I,J)=0.0 

8330  = 

DO  110  1=1, N 

8340= 

DO  no  J  =  1,M 

8350= 

110 

SYSM( I ,N+J)=-B( I , J) 

83G0=C 
8370  = 

DO  120  1=1, L 

8380= 

DO  120  J=1,N 

8390  = 
8400=C 

120 

SY3M(N+I, J)=C(I,J) 

8410  = 

DO  130  1=1, L 

8420= 

DO  130  J=1,M 

8430= 

130 

SYSM(N+I ,N+J)=D( I , J) 

8440=C 
8450=C 
8480  = 

DO  200  1=1, N 

8470  = 

DO  200  J=1,N 

to  00 
o  o 

II  II 

n 

200 

SKI,  J)=0.0 

8500  = 

DO  210  1=1, N 

8510  = 

210 

SI(I,I)=Z 

8520=C 

8530= 

CALL  NSUB(SI ,A,SIA,N,N) 

8540=C 

8550= 

DO  220  I=1,N 

8580= 

DO  220  J=1,N 

8570= 

8580=C 

220 

SYSMd,  J)=SIA(I,  J) 

8590= 

CALL  LS00F(SYSN,  12,L0,I10,U,  12 , 0 , 9 , MK  ,  lER ) 

8G00=C 

8610= 

8620=C 

IF(LO.GT.MO)  GO  TO  320 

8630  = 

DO  300  I=1,L0 

8640= 

IF<S( I ) .GT.EXl )  GO  TO  300 

8650= 

IF<S< I ) .LT.-EXl )  GO  TO  300 

8660= 

IF(S(I) .LT.EX2.AND.S( I ) .GT.-EX2)  IZ=1 

8670  = 
8680=C 

300 

CONTINUE 

8690= 

GO  TO  330 

8700=C 

8710= 

320 

DO  310  1=1, MO 

8720= 

IF(S( I ) .GT.EXl )  GO  TO  310 

8730= 

IF(S( I ) .LT.-EXl )  GO  TO  310 
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8740=  IF(S( I ) .LT.EX2.AND.S< I ) .GT.-EX2)  IZ=1 

8750=  310  CONTINUE 

8760=C 

70=  330  CONTINUE 

8780=C 

8780=  RETURN 

8800=  END 

8810=C 
8820=C 

8830=  SUBROUTINE  COMPLEX ( NN , ZERI » IFL) 

8840=  REAL  ZERKIB) 

8850=  INTEGER  IFLdS) 

8860=C 

8870=C  THIS  SUBROTINE  TESTS  FOR  COMPLEX  CONJIGATE  PAIRS. 

8880=C  A  COMPLEX  SYSTEM  ZERO  MILL  ALMAYS  OCCUR  IN  CONJIGATE  PAIRS 

8890=C  ANY  COMPLEX  ZERO  CALCULATED  WITHOUT  A  CONJIGATE  PAIR  WILL  BE 

B900=C  REJECTED  AS  A  SYSTEM  ZERO 

8910=C 

8920=  EX1=1.0E-S 

8930=C 

8940=  DO  100  1=1, NN 

8950=  IF( IFL( I ) .EQ.O)  GO  TO  100 

8960=  IF(2ERI(I).LT.EX1.AND.ZERI(I).GT.-EX1)  GO  TO  100 

8970=C 

8980=  DO  no  J  =  1,NN 

8990=  IF( IFL( J) .EQ.O)  GO  TO  110 

9000=  IF(ZERI( J).LT.EX1.AND.2ERI(J>.GT.-EX1>  GO  TO  110 

9010=C 

^0=  CCP=ZERI<I) 

^30=  CCN  =  ZERI<J) 

9040=  CCC=CCP+CCN 

9050=  IF(CCC.LT.EX1 .AND.CCC.GT.-EXl)  GO  TO  120 

9080=  GO  TO  130 

9070=  120  IFL(I)=1 

9080=  IFL(J)=1 

9090=  GO  TO  100 

8100=  130  IFL(I)=0 

9110=C 

9120=  no  CONTINUE 
9130=  100  CONTINUE 

9140=C 

9150=  RETURN 

9160=  END 

9170=*E0R 
9180**E0F 


Appendix  C 


Medn  Copg)Uter  Program 

This  appendix  contains  a  listing  of  the  computer  program  used  to 
design  controllers  and  simulate  them.  The  program  is  adopted  from  the 
program  used  by  Lt  Brett  Ridgely  (Ref  10) .  A  simplified  flow  chart  is 
shown  in  Figure  C-1.  The  program  listed  is  for  designs  with  actuator 
dynamics.  For  designs  without  actuator  actuator  dynamics  lines  2840 
through  3450  are  replaced  by  the  lines  3300  through  4380  shown  at  the 
end  of  the  main  program  listing. 

The  example  is  design  2A  simulated  at  the  design  condition  (5,000  ft, 
Mach  0.6)  and  also  at  the  condition  of  design  3  (10,000  ft,  Mach  0.9). 
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WELCOME  TO  THE  SINGULARLY  PERTURBED  ASYMPTOTIC  METHOD 
(SPAM  FOR  SHORT)  FOR  IRREGULAR  PLANTS  -  VERSION  3.2 
MAKE  SURE  BOTH  THE  IMSL  AND  CCBGOO  LIBRARIES  ARE  SET 

DESIGN  1  <1)  OTHER  (2)  .* 

450008  CM  STORAGE  USED 

1.812  CP  SECONDS  COMPILATION  TIME 
NON-FATAL  LOADER  ERRORS  - 
NON-EXISTENT  LIBRARY  GIVEN  -  SYS  10  1 

THE  All  MATRIX  IS  - 

0.0000 

THE  A12  MATRIX  IS  - 

0.0000  1.0000  .0318 

THE  A21  MATRIX  IS  - 

.0488 

0.0000 

0.0000 

THE  A22  MATRIX  IS  - 

-.3744  .0318  -.9SS5 

-39.9283  -2.1138  .0797 

S.2378  -.0063  -.7074 

THE  B1  MATRIX  IS  A  1  BY  3  ZERO  MATRIX 

THE  B2  MATRIX  IS  - 

.0531  .0046  .0250 

9.7929  -55.7388  2.6033 

-5.0544  -1.8927  4.3880 

SYSTEM  DIMENSION  =  4 
#  INPUTS/OUTPUTS  =  3 


ENTER  Cl  (  3  BY  1  ) 

ROW  1  :  0 
ROW  2  :  1 
ROW  3  :  0 

ENTER  C2  (  3  BY  3  > 

ROW  1  :  1,0,0 
ROW  2  :  0,0,0 
ROW  3  :  0,0,1 
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THIS  METHOD  IS  FOR  IRREGULAR  PLANTS,  THAT  IS,  FOR 
C2B2  NOT  OF  FULL  RANK 

***  TERMINAL  ERROR  (lER  =  129)  FROM  IMSL  ROUTINE  LUDATF 

TERMINAL  ERROR  (lER  =  130)  FROM  IMSL  ROUTINE  LIN03P 

THE  PLANT  IS  IRREGULAR 


ENTER  THE  M  MATRIX  (  3  BY  1  )  : 

ROW  1  :  0 
ROW  2  :  .25 
ROW  3  :  0 


THE  MATRIX  FI  IS  - 

0.0000 
1 . 0000 
0 . 0000 

THE  MATRIX  F2  IS  - 


1 . 0000  0 . 0000  0 . 0000 

0.0000  .2500  .0080 

0.0000  0.0000  1.0000 

THE  MATRIX  F2*B2  IS  - 

.0531  .0046  .0250 

2.4080  -13.9498  .5857 

-5.0544  -1.8927  4.3880 


NOW  VALUES  OF  SIGMA  MUST  BE  CHOSEN.  HOW  DO  YOU  WANT 
TO  INPUT  THESE  (0=DIAG0NAL,  1=FULL)  :  0 

ENTER  THE  3  VALUES  OF  SIGMA  (MUST  BE  POSITIVE)  I  1.15,30 
THE  MATRIX  KO  IS  - 

11.6944  .1981  -2.0615 

2.7389  -1.0522  -.1395 

14.6517  -.2256  4.4020 

THE  MATRIX  K1  MUST  NOW  BE  INPUT.  HOW  DO  YOU  WISH 
TO  DO  THIS  (0=SCRATCH  ,  1=MULTIPLE  OF  KO)  :  1 

WHAT  MULTIPLE  OF  KO  IS  K1  I  2 

THE  K1  MATRIX  IS  - 

23.3888  .3962  -4.1231 

5.4778  -2.1044  -.2789 

29.3035  -.4513  8.3041 


117 


HE  FIRST  SET  OF  SLOW  EIGENVALUES  ARE 


LAMBDA(l)  =  (-2..0.) 
LAMBDA (2)  =  (-2.,0.) 


LAMBDAO)  =  (-2.»0.) 


THE  SECOND  SET  OF  SLOW  EIGENVALUES  ARE  I 
LAMBDA(4)  =  (-4.,0.) 

THE  FAST  EIGENVALUES  ARE  (MUST  BE  MULTIPLIED  BY  G)  : 
LAMBDA(5)  =  (-l.,0.) 


LAMBDA (6)  =  (-30.,0.) 

LAMBDA(7)  =  (-15., 0.) 

FORM  CLOSED  LOOP  MATRICES  (1),  CHANGE  DUTPUTS/EIGENVALUES 
(2) ,  OR  STOP  (0)  :  1 

PRINT  C-L  MATRICES  (1=YES)  I  0 

ENTER  THE  GAIN  CONSTANT  ,  G  :  29 


ENTER  THE  3  ACTUATOR  CONTANTS  : 

RUDDER  CONSTANT  .*  20 
AILERON  CONSTANT  :  20 
CANNARD  CONSTANT  I  20 


THE  OVERALL  CLOSED  LOOP  EIGENVALUES  ARE  I 


EIGENVALUE 

(  1  ) 

EIGENVALUE 

(2) 

EIGENVALUE 

(3) 

EIGENVALUE 

(4) 

EIGENVALUE 

(5) 

EIGENVALUE 

(6) 

EIGENVALUE 

(7) 

EIGENVALUE 

(8) 

EIGENVALUE 

(9) 

EIGENVALUE 

(  10) 

(-9.34B964B74272, 131 .512SG0S02B) 
(-9.348964874273,-131 .512SG0B02B) 
(-8.041546031385,92.59541875491 ) 
(-8.041546031385,-92.59541875491 ) 
(-9.135660521417,21.57045712481 ) 
(-9.135660521417,-21.57045712481 ) 
(-4.027317394378,0. ) 
(-2.104785028134,0.  ) 
(-2.010762726775,0. ) 

(-2.000431356545,0. ) 


CHANGE  OUTPUTS/EIGENVALUES  (1),  CHANGE  ONLY  GAIN  (2), 
CHANGE  SIGMA  (KO/Kl)  (3),  INTEGRATE  C-L  EONS  (4), 
CHANGE  PLANT  (5>,  STOP  (0)  14 

*•»•»■»■»■#  ENTER  INTEGRATING  ROUTINE 
FOR  THIS  SYSTEM  THERE  ARE  3  INPUTS  AND  3  OUTPUTS 
DO  YOU  WISH  TO  SEE  THE  INPUT  HISTORIES  OR  OUTPUT 
HISTORIES  ?  ( 1=INPUT,2=0UTPUT)  I  2 
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ENTER  MANEUVER  ( Y  .  P .  =  1  ,  W .  L .  T .  =2 ,  H  .  T  .  =3  )  .*  1 
ENTER  INPUT  COMMAND 

FINAL  VALUE,  RAMP  TIME  (UNIT  STEP  =  1,0) 

1  :  -3,1 

2  :  0,0 


ENTER  STEP  RESPONSE  END 
3  :  3,0 

ENTER  STEP  RESPONSE  END 

ENTER  FINAL  TIME  AND  TIME 

T  OUTPUTS 

time  (INFINITY 

TIME  (INFINITY 

STEP  :  1.75,. 05 

o 

■  t  ■  • 

o  o 

II  11 

0 . 00 

0 . 0000 

0.0000 

0 . 0000 

0 . 0000 

0 . 0000 

.05 

-.1936 

.0014 

1  .  1608 

.  0580 

-.1356 

.  10 

-.3831 

.  0028 

1.9707 

.  1566 

-.2256 

.  15 

-.5369 

.  0044 

2.4866 

.2809 

-.2560 

.20 

-.6728 

.0062 

2.7875 

.4203 

«  w'  w 

.25 

-.8110 

.0073 

2.9457 

.5676 

-.2434 

.30 

-.9563 

.  0095 

3.01 73 

.713^ 

-.2379 

.35 

-1 . 1077 

.0112 

3.0411 

.  8  705 

-.2372 

.40 

-1.2591 

.0129 

3.0414 

1.0226 

-.2365 

.  45 

-1 .4072 

.0144 

3.0323 

1 .1742 

-.2330 

.  50 

-1.5562 

.0159 

3.0213 

1.3252 

-.2310 

.55 

-1.7034 

.0174 

3.0115 

1.4758 

-.2276 

.60 

-1.8586 

.0188 

3.0040 

1.6260 

-.2326 

.65 

-2.0098 

.0201 

2.9999 

1.7759 

-.2339 

.70 

-2.1550 

.0212 

2.9955 

1.9257 

-.2293 

.75 

-2.3016 

.0223 

2 . 9934 

2.0754 

-.2262 

.80 

-2.4463 

.0233 

2.9923 

2.2250 

-.2213 

.85 

-2 . 6002 

.0242 

2.9916 

2.37*6 

- . 2256 

.SO 

-2.7431 

.0251 

7  aci  i  7 

2 . 5242 

-.22*9 

.95 

-2. 3376 

.0259 

2.9910 

2.6737 

-.2239 

1 . 00 

-3.0552 

.0267 

2.9909 

2.8233 

-.2319 

1.05 

-3.0155 

.0261 

1.7695 

2.9117 

- . 1038 

1  .  10 

-2.9676 

.0254 

1.2151 

2 . 9725 

.0049 

1  . 15 

-2.9602 

.02^1 

,7618 

3.0106 

.0504 

1.20 

-2.9779 

.0229 

.4358 

3.0324 

.0545 

1.25 

-2.9932 

.0217 

.2229 

3.0435 

.  0503 

1.30 

-2.9964 

.0205 

.0954 

3.0483 

.0519 

1.35 

-2.9930 

.0192 

.0261 

3.0496 

.0566 

1.40 

-2.9903 

.  0 1 78 

-.0069 

3.0492 

.0589 

1 . 45 

-2.9907 

.0166 

-.0193 

3.0483 

.0576 

1.50 

-2.9925 

.0153 

-.0211 

3.0472 

.0547 

1.55 

-2.9940 

.0142 

-.0183 

3.0463 

.0523 

1.60 

-2.9946 

.0131 

-  .01'i2 

3.0456 

.0510 

1.65 

-2 . 99dg 

.0120 

-.010* 

3.0451 

.  0502 

1.70 

-2.9352 

.011 0 

- .0073 

3.0447 

.0495 

1.75 

-2.9357 

.0101 

- . 005 1 

3 . 0444 

.  0488 

BETA  COMMAND  =-1.6330318^6067 
R  COMMAND  =2.823251138512 
PHI  =.02907321538383 
BETA  TRANS.  =.011736696361^2 
R  TRANS.  =.3810131352851 
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HANGE  OUTPUTS/EIGENi^ALUES  <1),  CHANGE  ONLY  GAIN  <2), 

CHANGE  SIGMA  (KO/Kl)  (3),  INTEGRATE  C-L  EQNS  (4), 

CHANGE  PLANT  (5),  STOP  (0)  :  5 

NEW  DESIGN  (1)  OR  ROBUSTNESS  TEST  (2)  ?!  2 

ENTER  THE  DIMENSION  OF  THE  SYSTEM  :  4 
ENTER  THE  NUMBER  OF  ACTIVE  INPUTS  (SIZE  OF  B2)  :  3 
ENTER  All  (  1  BY  1  ) 

ROW  1  :  0 

ENTER  A12  (  1  BY  3  ) 

ROW  1  :  0,1,. 0255 

ENTER  A21  (  3  BY  1  ) 

ROW  1  :  .0332 
ROW  2  :  0 
ROW  3  :  0 

ENTER  A22  (  3  BY  3  ) 

ROW  1  :  -.4836, .0255,-. S9S7 
ROW  2  :  -75.3125,-3.9173, .0382 
ROW  3  :  10. 6995,-. 0294, -.5117 

Bl  IS  A  ZERO  MATRIX  OF  DIMENSION  1  BY  3 

ENTER  B2  (  3  BY  3  ) 

ROW  1  :  .0473, .0066, .0399 

ROW  2  :  13.4743,-82.0763,8.3667 

ROW  3  :  -7.7774,-3.2921,8.7542 

***  TERMINAL  ERROR  ( lER  =  129)  FROM  IMSL  ROUTINE  LUDATF 

*■»*  TERMINAL  ERROR  (lER  =  130)  FROM  IMSL  ROUTINE  LINV3F 

THE  PLANT  IS  IRREGULAR 


THE  MATRIX  FI  IS  - 

0 . 0000 
1 . 0000 
0.0000 

THE  MATRIX  F2  IS  - 

1.0000  0.0000  '0.0000 

0 . 0000  . 2500  . 006  4 

0 . 0000  0 . 0000  1 . 0000 
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HE  MATRIX  F2-»B2  IS 


f? 


.0473  .0066  .03SS 

3.3190  -20.5401  2.1530 

-7.7774  -3.2921  8.7542 

THE  MATRIX  KO  IS  - 

11.6944  .1981  -2.0615 

2.7389  -1.0522  -.1395 

14.6517  -.2256  4.4020 

THE  K1  MATRIX  IS  - 

23.3888  .3962  -4.1231 

5.4778  -2.1044  -.2799 

29.3035  -.4513  8.8041 

THE  FIRST  SET  OF  SLOW  EIGENVALUES  ARE  ! 

LAMBDA(l)  =  (-2.r0.) 

LAMBDA(2)  =  (-2.,0.) 

LAMBDA(3)  =  (-2.,0.) 

THE  SECOND  SET  OF  SLOW  EIGENVALUES  ARE  : 

LAMBDA (4)  =  <-4.,0.) 

THE  FAST  EIGENVALUES  ARE  (MUST  BE  MULTIPLIED  BY  G) 
LAMBDA(5)  =  (-1.11882838742,0.) 

LAMBDA(6)  =  (-55.06024604*91,0.) 

LAMBDA(7)  =  (-21.78931709532,0.) 


THE  OVERALL  CLOSED  LOOP  EIGENVALUES  ARE  : 


EIGENVALUE  (1) 
EIGENVALUE  (2) 
EIGENVALUE  (3) 
EIGENVALUE  (4) 
E  - GE^^iVALUE  (5) 
EIGENVALUE  (6) 
EIGENVALUE  (7) 
EIGENVALUE  (S) 
E I G  E  N  ^r‘  A  L  U  E  (  9  ) 
EIGENVALUE  (  iO 


=  (-9.511823000371 , 176.4317591339) 

=  (-9.511823000371,-178.4317581339) 
=  (-8.918649595212,1 1 1.9021929101 ) 


<'-8.91 96495952 12,-111. 902  1 929 1 0  ■  ) 
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NTER  INTEGRATING  ROUTINE  ****** 

FOR  THIS  SYSTEM  THERE  ARE  3  INPUTS  AND  3  OUTPUTS 
DO  YOU  WISH  TO  SEE  THE  INPUT  HISTORIES  OR  OUTPUT 
HISTORIES  ?  ( 1=INPUT,2=0UTPUT)  :  2 
ENTER  MANEUVER  ( Y . P . = 1 , W . L . T . =2 , H . T . =3 )  :  1 


ENTER  INPUT  COMMAND 

FINAL  VALUE,  RAMP  TIME  (UNIT  STEP  ^  1,0) 

1  :  -3,. 75 

2  :  0,0 

ENTER  STEP  RESPONSE  END  TIME  (INFINITY=0)  :  0 
3  :  4,0 

ENTER  STEP  RESPONSE  END  TIME  (INFINITY=0)  :  .75 
ENTER  FINAL  TIME  AND  TIME  STEP  :  1 . 5 , . 05 


T  OUTPUTS 


0.00 

0.0000 

0.0000 

0 . 0000 

0 . 0000 

0 . 0000 

.05 

-.2520 

-.0216 

6.0860 

.3043 

.0523 

.  10 

-.5291 

-.0078 

3.2783 

.4682 

-.0609 

.  15 

-.7149 

-.0059 

4.0252 

.6695 

-.0454 

.20 

-.8935 

.0076 

4.2972 

.8843 

- . 009 1 

.25 

-1.0717 

.0090 

3.6887 

1.0688 

-.0029 

.30 

-1.2708 

.0118 

4.2317 

1.2804 

.0096 

.35 

-1.4809 

.0139 

3.8643 

1.4736 

-.0074 

.40 

-1.6837 

.0157 

4.0574 

1.6764 

-.0073 

.45 

-1.8823 

.0184 

3.9882 

1.8759 

- . 0065 

.50 

-2.0740 

.0201 

3 . 9926 

2.0755 

.0015 

.55 

-2.2712 

.0218 

4.0049 

2.2757 

.  0045 

.60 

-2.4731 

.0234 

3.9821 

2.4748 

.0017 

.65 

-2.6779 

.0253 

3.9960 

2.6746 

-.0033 

.70 

-2.8755 

.0273 

3.9971 

2.8745 

-.0010 

.75 

-3.0723 

.0297 

3.9862 

3.0738 

.0015 

.80 

-3.0254 

,0543 

-1 .8588 

2.9809 

-.0443 

.85 

-2.9424 

.0390 

.4810 

3.0049 

.0625 

.90 

-2.9556 

.0379 

.  1216 

3.01 10 

.  055A 

.95 

-2.9753 

.  0255 

-.3879 

2.9916 

.01S3 

1.00 

-2.9935 

.0247 

.3308 

3.0031 

.0147 

1  .05 

-2.9927 

.0231 

-.2394 

2.9962 

.0034 

1  .  10 

-2.9859 

.0216 

.  1051 

3.0014 

.0155 

1  .  15 

-2.9858 

.0209 

-.0525 

2.9938 

.  0 1 30 

1.20 

-2.9877 

.  0 1 90 

-.0132 

2.9981 

.010^ 

1 .25 

-2.9911 

.0180 

.0088 

2.9986 

.0075 

1.30 

-2.9922 

.0166 

,0263 

2.9973 

.  0050 

1.35 

-2.9925 

.0154 

.0059 

2.9976 

.0050 

1.40 

-2.9929 

.0141 

-.0141 

2.9968 

.0039 

1.45 

-2.9936 

.0128 

- .  0033 

2.9967 

.0031 

1.50 

-2.9944 

.0118 

-.0058 

2.9964 

.0019 

BETA  COMMAND  =-1.261148520343 
R  COMMAND  =3.073795857869 
PHI  =.03064685696701 
BETA  TRANS.  =.01610555728732 
R  TRANS.  =.3820396586508 
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100  = 
110= 
V20= 
J0  = 
140= 
150= 
1B0= 
170= 
180= 
190= 
200= 
210  = 
220= 
230= 
240= 
250= 
2B0= 
270= 
280= 
290= 
300  = 
310= 
320= 
330  = 
340  = 
350  = 
3B0= 
370  = 
^0  = 
\xd0  = 
400  = 
410= 
420= 
430= 
440  = 
450= 
4B0  = 
470= 
480  = 
490  = 
500= 
510= 
520= 
530  = 
540  = 
550= 
580  = 
570  = 
580  = 
590  = 
800  = 
810  = 
820  = 
F30= 


PROGRAM  PORTER( INPUT, OUTPUT, TAPE5=INPUT,TAPEB=0UTPUT, 

1  PL0T,TAPE4=PL0T) 

C 

REAL  All ( 10, 10) ,A12< 10, 10) ,A21 (10,10) ,AZ2( 10, 10) 

REAL  B1 ( 10, 10) ,B2( 10, 10) ,C1 ( 10, 10) ,C2( 10, 10) 

REAL  FI ( 10, 10) ,F2( 10, 10) ,M(10,10) ,UC(10) 

REAL  K0(10,10),K1(10,10),X(10),Y(10),Z(10) 

REAL  A12FF( 10, 10) ,0< 10) ,C2B2( 10, 10) 

REAL  WKAREA( 10) ,MA1 1(10,10) , MA12 ( 10 , 10 ) , F2B2 ( 10 , 10 ) 

REAL  SIGMA(10,10) ,FBINU( 10, 10) ,K00(10,10) ,K11(10,10) 

REAL  BETA ( 10) , 3L0W2 ( 10 ) , AFIF ( 10 , 10 ) , F2IF1 ( 10 , 10 ) , F2IN0 ( 10 , 10 ) 
COMPLEX  ALFA( 10) , LAMBDA ( 10) ,LAM( 10) , W( 10) , W1 ( 10) ,ZQ( 10, 10) 

REAL  AY ( 10, 10) , BEE ( 10, 10) , CEE ( 10, 10) 

REAL  AC ( 10, 10) ,ACK1 ( 10,10) ,ACKF1 (10,10) , ACKF2 ( 10 , 10 ) ,ACKO( 10, 10) 
REAL  FB2( 10, 10) ,F2B2K0( 10, 10) ,FAST( 10, 10) 

REAL  IDEN( 10, 10) ,F2S( 10, 10) ,AY1 (10,10) 

REAL  XP( 10) ,MULT,ZDaT( 10) ,U( 10) ,ER(5) ,RR(2) 

REAL  VC( 10) ,ODT( 10) , SLOPE ( 10) ,OST( 10) 

INTEGER  KOR(  10)  ,K'.»ST(  10) 

DIMENSION  UI0RK(500)  ,IMORK(10) 

COMMON  AY,BEE,0,NL,L 
EXTERNAL  F 
C 
C 

DO  99  1=1,3 
99  PRINT*,"  " 

C 

PRINT*, "WELCOME  TO  THE  SINGULARLY  PERTURBED  ASYMPTOTIC  METHOD" 
PRINT*, "(SPAM  FOR  SHORT)  FOR  IRREGULAR  PLANTS  -  OERSION  3.2" 
PRINT*, "MAKE  SURE  BOTH  THE  IMSL  AND  CCBBOO  LIBRARIES  ARE  SET" 
PRINT*,"  " 

C 

X10=0 
Xll  =  l 

PRINT*,"  DESIGN  1  (1)  OTHER  (2) 

READ* , XX 

IF(XX.Ea.l)  GOTO  10 
C 

9970  IF(X10.NE.5)  GO  TO  9971 
PRINT*,"  " 

PRINT*,"  NEW  DESIGN  (1)  OR  ROBUSTNESS  TEST  (2)  ?  I" 

READ*,X11 

9971  CONTINUE 
PRINT*,"  " 

PRINT*, "ENTER  THE  DIMENSION  OF  THE  SYSTEM  !  " 

READ*,N 

C 

IF  (N.LE.7)  GOTO  100 

PRINT*, "THE  DIMENSIONS  MUST  BE  CHANGED  TO  HANDLE  THAT" 

GOTO  9999 
C 

100  PRINT*, "ENTER  THE  NUMBER  OF  ACTIOE  INPUTS  (SIZE  OF  B2)  :  " 
READ*,L 
NM=N-L 
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640=C 
650= 
660= 
-70  = 
:,80  = 
690  = 
700= 
710  = 
720= 
730= 
740= 
750= 
760=C 
770= 
780= 
790=C 
800= 
810= 
820  = 
830=C 
840= 
850  = 
860  = 
870=6 
880= 
890=C 
900= 


940  = 
950= 
960= 
970= 
980  = 
990=C 
1000= 
1010= 
1020= 
1030= 
1040»C 
1050= 
1060= 
1070= 


200 


10 


30 


20 


1080= 

1090= 

1100= 

1110=C 

1120= 

1130= 

1140= 

1150=  9972 
1160= 

1170= 


PRINT* r"  ENTER  All  ( 
CALL  READM(NM,NM,A11) 
PRINT*."  " 

PRINT*,"  ENTER  A12  ( 
CALL  READM<NI1,L,A12) 
PRINT*,"  " 

PRINT*,"  ENTER  A21  ( 
CALL  REA0M(L,Nri,A21  ) 
PRINT*,"  " 

PRINT*,"  ENTER  A22  ( 
CALL  READM(L,L,A22) 


",NM,"  BY  ",NM,"  )" 

",NM,"  BY  ",L,"  )" 

",L,"  BY  ",NM,"  )" 

",L,"  BY  ",L,"  )" 


PRINT*,"  " 

PRINT*, "B1  IS  A  ZERO  MATRIX  OF  DIMENSION  ",NM,"  BY  ",L 


DO  200  I=1,NM 
DO  200  J=1,L 
Bid,  J)=0. 


PRINT*,"  " 

PRINT*,"  ENTER  B2  (  ",L,"  BY  ",L,"  )" 
CALL  READM<L,L,B2) 

IF(X11.EQ.2)  GO  TO  9972 


GOTO  20 


CALL  DSGNl <N,L,A1 1 ,A12,A21 ,A22,B1 ,B2) 
PRINT*,"  " 

PRINT*,"  SYSTEM  DIMENSION  =  ",N 
PRINT*,"  #  INPUTS/OUTPUTS  =  ",L 
PRINT*,"  " 

NM=N-L 
GOTO  20 


PRINT*,"  " 

PRINT*,  "NEW  OUTPUT  MATRICES  ?<1  =  YES) 

READ*,X5 

IF  (X5.NE.1)  GOTO  380 
PRINT*,"  " 

PRINT*,"  ENTER  Cl  (  ",L,"  BY  ",NM,"  )" 

CALL  READM(L,NM,C1 ) 

PRINT#,"  " 

PRINT*,"  ENTER  C2  (  ",L,"  BY  ",L,"  )" 

CALL  READM(L,L,C2) 

PRINT*,"  " 

PRINT*, "THIS  METHOD  IS  FOR  IRREGULAR  PLANTS,  THAT  IS,  FOR" 
PRINT*, "C2B2  NOT  OF  FULL  RANK" 

CALL  UMULFF(C2,B2,L,L,L, 10, 10,C2B2, 10, lER) 

Dl  =  l  . 

CALL  LIN03F(C2B2,B,4,L,10,D1,D2,WKAREA, lER) 
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1180*  DET=D1*2**D2 

1190=  IF  (DET.EQ.O)  GOTO  280 

1200=  PRINT*, "THIS  PLANT  IS  REGULAR  -  PROGRAM  TERMINATED" 

,V^10=  GOTO  9999 

■-  .20=  280  PRINT#,  "THE  PLANT  IS  IRREGULAR" 

1230=  PRINT*,"  " 

1240=  IF(X11.EQ.2)  GO  TO  9973 

1250=C 

12B0=  380  CONTINUE 

1270=  810  PRINT*,"  " 

1280=  PRINT*,"  " 

1290=  PRINT*, "ENTER  THE  M  MATRIX  (  ",L,"  BY  ",NM,"  )  I  " 

1300=  CALL  READM(L,NM,M) 

1310=  PRINT*,"  " 

1320=  9973  CALL  VMULFF ( M , A1 1 , L , NM , NM , 10 , 10 , MAI  1 , 10 , lER ) 

1330=  CALL  MADD(C1,MA11,L,NM,F1) 

1340=  CALL  yMULFF(M,A12,L,NM,L, 10, 10,MA12, 10, lER) 

1350=  CALL  MADD(C2,MA12,L,L,F2> 

1380=  CALL  MADD(C2,MA12,L,L,F2S) 

1370=  CALL  0MULFF<F2,B2,L,L,L,10,10,F2B2,10,IER) 

1380=  PRINT*,"  " 

1390=  PRINT*, "THE  MATRIX  FI  IS  -" 

1400=  CALL  PRINTM(L,NM,F1 ) 

1410=  PRINT*,"  " 

1420=  PRINT*, "THE  MATRIX  F2  IS 

1430=  CALL  PRINTM(L,L,F2) 

1440=  PRINT*,"  " 

1450=  PRINT*, "THE  MATRIX  F2*B2  IS  -" 

14G0»  CALL  PRINTM(L,L,F2B2) 

^0=  Dl  =  l 

1^0=  CALL  LIN03F(F2B2,B,4,L,10,D1,D2,WKAREA,IER) 

1490=  DET=D1*2**D2 

1500=  IF  (DET.NE.O)  GOTO  370 

1510=  PRINT*, "THE  MATRIX  F2*B2  IS  NOT  FULL  RANK  -  TRY  AGAIN" 

1520=  GOTO  30 

1530=  370  CONTINUE 

1540=  IF(Xll.Ea.2)  GO  TO  9974 

1550=  PRINT*,"  " 

1560=  PRINT*, "NUW  VALUES  OF  SIGMA  MUST  BE  CHOSEN.  HOW  DO  YOU  WANT" 

1570=  PRINT*, "TO  INPUT  THESE  (0=DIAG0NAL,  1=FULL)  I" 

1580=  READ*,JJ 

1590=  IF  <JJ.EQ.1>  GOTO  390 

1G00=C 

1S10=  DO  393  1=1, L 

1620=  DO  393  J=1,L 

1630=  393  SIGMA(I,J)=0 

1640=C 

1650=  PRINT*,"  " 

1680=  PRINT*,  "ENTER  THE  ",L,"  VALUES  OF  SIGMA  (MUST  BE  POSITIVE) 

1B70=C 

1680=  DO  400  I=1,L 

1690=  400  REAO*,SIGMA<I,I) 

1700»C 

1710=  GOTO  410 
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1720=C 

1730*  390  PRINT#,"  " 

1740*  PRINT#, "ENTER  SIGMA 

*^.  *50*  CALL  READM(L,L, SIGMA) 

i :o0*C 

1770=  9974  CONTINUE 

1790*  410  CALL  0MULFF(F2,B2,L,L,L, 10, 10,F2B2, 10, lER) 

1790*C 

1800*  DO  392  1*1, L 

1810*  DO  392  J*1,L 

1820*  392  FB2<I, J)*F2B2<I, J) 

1830*C 

1840*  IF(X11.EQ.2)  GO  TO  9975 

1850=  CALL  LIN01F(F2B2,L, 10,FBINO, 10,WKAREA, lER) 

1860*  CALL  OMULFF(FBINO,SIGMA,L,L,L, 10, 10,K0, 10, lER) 

1870=C 

1880=  9975  PRINT#,"  " 

1890*  PRINT#, "THE  MATRIX  KO  IS 

1900*  CALL  PRINTM<L,L,KO) 

1910*C 

1920*  DO  4B0  I=1,L 

1930*  DO  460  J*1,L 

1940*  460  K00( I , J)=-K0< I , J) 

1950*C 

1960*  IF(Xll.Ea.2)  GO  TO  9976 

1970*  PRINT#,"  " 

1980*  PRINT#, "THE  MATRIX  K1  MUST  NOW  BE  INPUT.  HOW  DO  YOU  WISH" 

1990*  PRINT#, "TO  DO  THIS  (0*SCRATCH  ,  1=MULTIPLE  OF  KO)  I" 

2000*  READ#,G2 

Vo*  9976  IF  (G2.Ea.l)  GO  TO  480 

2020*  520  CONTINUE 

2030*  PRINT#,"  " 

2040*  PRINT#,"  ENTER  K1  <",L,"  BY  ",L,")" 

2050*  CALL  READM(L,L,K1 ) 

2060*C 

2070*  GOTO  510 

2080*  480  IF  <X11.EQ.2)  CO  TO  9977 

2090*  PRINT#,"  " 

2100*  PRINT#, "WHAT  MULTIPLE  OF  KO  IS  K1  I" 

2110*  READ#, MULT 

2120*C 

2130*  DO  522  1*1, L 

2140*  DO  522  J*1,L 

2150*  522  Kid,  J)=MULT#KO(I,J) 

2160*C 

2170*  9977  PRINT#,"  " 

2180*  PRINT#,"  THE  K1  MATRIX  IS  -" 

2190=  CALL  PRINTM(L,L,K1  ) 

2200*C 

2210*  DO  523  1*1, L 

2220*  DO  523  J*1,L 

2230*  523  K11(I,J)*K1(I,J) 

2240*C 

2250-  510  CALL  EIGZF(K 1 1 , 10 , KOO , 10 ,L,0, ALFA, BETA , ZQ , 10 , WKAREA , lER ) 
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22S0=C 

2270=  DO  530  I=1,L 

2280=  YY=BETA(I) 

r:  '30=  IF  (YY.EQ.O)  GOTO  540 

2cb0=  LAMBDA(I)=ALFA<I)/BETA(I) 

2310=  GOTO  530 

2320=  540  LAMB0A(I)=9.E+99 

2330=  530  CONTINUE 

2340=C 

2350=  PRINT*,"  " 

2360=  PRINT*, "THE  FIRST  SET  OF  SLOW  EIGENVALUES  ARE  I" 

2370=  PRINT*,"  " 

2380=C 

2390=  DO  550  I=1,L 

2400=  PRINT*, "LAMBDA( " , I , " )  =  ",LAMBDA(I) 

2410=  550  PRINT*,"  " 

2420=C 

2430=  CALL  LINVIF ( F2S ,L , 10 , F2INV, 10 , WKAREA , lER ) 

2440=  CALL  VMULFF(F2INV,F1 ,L,L,NM, 10, 10,F2IF1 , 10, lER) 

2450=  CALL  VMULFF( A12,F21F1 ,NM,L,NM, 10, 10, A12FF, 10, lER) 

24S0=C 

2470=  DO  560  I=1,NM 

2480=  DO  560  J=1,NM 

2490=  560  AFIF<I,J)=-A12FF(I,J) 

2500=C 

2510=  CALL  MADD(A11 ,AFIF,NM,NM,SL0W2) 

2520=  CALL  EIGRF(SL0W2,NM, 10,0,M,ZQ, 10, WKAREA, lER) 

2530=C 

^0=  PRINT*,"  " 

A50=  PRINT*,  "THE  SECOND  SET  OF  SLOW  EIGENVALUES  ARE 

2560=  PRINT*,"  " 

2570=C 

2580=  DO  580  I=1,NM 

2590=  J=I+L 

2600=  PRINT*, "LAMBDA (",J,")  =  ",W(I) 

2610=  580  PRINT*,"  " 

2620=C 

2630=  CALL  VMULFF(FB2,K0,L,L,L, 10, 10,F2B2K0, 10, lER) 

2640=C 

2650=  DO  524  I=1,L 

2660=  DO  524  J=1,L 

2670=  524  FAST<I,J)=-F2B2K0(I,J) 

2680=C 

2690=  CALL  EIGRF(FAST,L, 10,0, W1 ,ZQ, 10, WKAREA, lER) 

2700=C 

2710=  PRINT*, "THE  FAST  EIGENVALUES  ARE  (MUST  BE  MULTIPLIED  BY  G)  I" 

2720=  PRINT*,"  " 

2730=  DO  590  I=1,L 

2740=  J=NM+L+I 

2750=  PRINT*, "LAMBDA (",J,")  =  ",W1(I) 

2760=  590  PRINT*,"  " 

2770=  IF(Xll.Ea.2)  GO  TO  9978 

2780=C 

2790=  PRINT*, "FORM  CLOSED  LOOP  MATRICES  (1),  CHANGE  OUTPUTS/EIGENVALUES" 
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2800=  PRINT*,"  (2),  OR  STOP  (0) 

2810=  READ«,X3 

2820=  IF  {X3.Ea.2)  GOTO  30 

2PgO=  IF  (X3.EQ.0)  GOTO  3399 

..  ;;0=  3999  PRINT#,"  " 

2850=  PRINT*, "PRINT  C-L  MATRICES  (1=YES)  I" 

28S0=  READ*,X25 

2870=  PRINT*,"  " 

2880=  PRINT*, "ENTER  THE  GAIN  CONSTANT  ,  G  I" 

2890=  REAO*,G 

2900*C 

2910=  DO  325  1  =  1, L 

2920=  DO  325  J=1,L 

2930=  325  AC(I,J)=0 

2940=C 

2950=  PRINT*,"  " 

2960=  PRINT*, "ENTER  THE  ",L,"  ACTUATOR  CONSTANTS  I" 

2970=  PRINT*, "RUDDER  CONSTANT  :" 

2980=  READ*, AC (1,1) 

2990=  PRINT*, "AILERON  CONSTANT  I" 

3000=  READ*, AC (2, 2) 

3010=  PRINT*, "CANNARD  CONSTANT  I" 

3020=  READ*, AC (3, 3) 

3030=C 

3040=  9978  CONTINUE 
3050=  I1=N+L+L 

3060=  I2=L+NM 

3070=  I3=L+L+NM 

3080=C 

«0=  DO  600  1  =  1,11 

3T00=  DO  620  J=1,I1 

3110=  620  AY(I,J)=0.0 

3120=  DO  630  K=1,L 

3130=  BEE(I,K)=0.0 

3140=  630  CEE(K,I)=0.0 

3150=  600  CONTINUE 

3160=C 

3170=  CALL  0MULFF(K0,F2,L,L,L,10,10,K0F2,10,IER) 

3180=  CALL  VMULFF<AC,K1 ,L,L,L, 10, 10,ACK1 , 10, lER) 

3190=  CALL  VMULFF(AC,K0F2,L,L,L, 10, 10,ACKF2, 10, lER) 

3200=  CALL  VMULFF<AC,K0,L,L,L,10,10,ACK0,10,IER) 

3210=  CALL  0MULFF(ACK0,F1 ,L,L,NM, 10, 10,ACKF1 , 10, lER) 

3220»C 

3230=  DO  640  I=1,NM 

3240=  DO  640  J=1,NM 

3250=  640  AY(L+I,L+J)=A11<I, J) 

3260=C 

3270=  DO  650  I=1,L 

3280=  DO  650  J=1,NM 

3290*  AY(I,L+J)*-F1<I,J) 

3300=  AY( I2+I ,L+J)»A21 ( I , J) 

3310*  AY<I3+I,L+J)»-G*ACKF1(I,J) 

3320=  AY(L+J,I2+I)»A12( J,I) 

3330=  650  CEE(I,L*J)=C1(I, J) 
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3340=C 

3350-  DO  660  1=1 rL 

3360=  DO  660  J=lfL 

r.''70»  AY(I,I2+J)=-F2(I,J) 

1-80=  AY( I2+I , I2+J)=A22( I r J) 

3390=  AY(I2+I,I3+J)=B2(I, J> 

3400=  AY(I3+I, J)=G*ACK1{I, J) 

3410=  AY< I3+I , I2+J ) =-G*ACKF2( I , J) 

3420=  AY(I3+I»I3+I)=-AC(I,I) 

3430=  CEE(I,I2+J)=C2(I, J) 

3440=  BEE<I»I)=1.0 

3450=  660  BEE(I3+Ir J)=G*ACKO<Ir J) 

3460=0 

3470=  NL=I1 

3480=  IF  (X25.NE.1)  GOTO  4310 

3490=  PRINT*."  " 

3500=  PRINT*. "THE  CL08ED  LOOP  AY  MATRIX  IS 

3510=  CALL  PRINTM(NL.NL,AY) 

3520=  PRINT*."  " 

3530=  PRINT*. "THE  CLOSED  LOOP  BEE  MATRIX  IS  -" 

3540=  CALL  PRINTM(NL.L.BEE) 

3550=  PRINT*."  " 

3560=  PRINT*. "THE  CLOSED  LOOP  CEE  MATRIX  IS 

3570=  CALL  PRINTM ( L . NL . CEE ) 

3580=  4310  CONTINUE 
3590=  PRINT*."  " 

3S00=C 

3610=  DO  4620  I=1,NL 

3R20=  DO  4620  J=1,NL 

<5*10=  4620  AY1(I.J)=AY(I.  J) 

3640=C 

3650=  PRINT*. "THE  OVERALL  CLOSED  LOOP  EIGENVALUES  ARE  I" 

3660=  PRINT*."  " 

3670=  CALL  EIGRF ( AYl . NL . 10 . 0 . HI . ZQ . 10 . MKAREA . lER ) 

3680=C 

3690=  DO  4700  I=1.NL 

3700=  4700  PRINT*."  EIGENVALUE  (".I.")  =  ".WKD 

3710=  IF(X11.EQ.2)  GO  TO  9979 

3720=C 

3730=  4710  CONTINUE 
3740=  PRINT*."  " 

3750=  PRINT*. "CHANGE  OUTPUTS/EIGENVALUES  (1).  CHANGE  ONLY  GAIN  (2)." 

3760=  PRINT*. "CHANGE  SIGMA  <K0/K1)  <3).  INTEGRATE  C-L  EQNS  (4)." 

3770=  PRINT*."  CHANGE  PLANT  <5).  STOP  (0)  :" 

3760=  READ*.X10 

3790=  IF  (XIO.EQ.O)  GOTO  9999 

3800=  IF  (XlO.EQ.l)  GOTO  30 

3810=  IF  (X10.EQ.2)  GOTO  3999 

3820=  IF  (X10.EQ.3)  GOTO  370 

3830=  IF  (XIO.EQ.S)  GOTO  9970 

3840=C 
3850=C 

3880=  9979  PRINT*."  " 

3870=  PRINT*."  ******  ENTER  INTEGRATING  ROUTINE  *****»" 
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3880=  PRINT*, "FOR  THIS  SYSTEM  THERE  ARE  ",L,"  INPUTS  AND  "rL,"  OUTPUTS" 

3890=  PRINT*, "DO  YOU  WISH  TO  SEE  THE  INPUT  HISTORIES  OR  OUTPUT" 

a900=  PRINT*,"  HISTORIES  ?  ( 1 = INPUT , 2=0UTPUT )  :" 

\\0-  READ*,AA1 

3920=  PRINT*,"  ENTER  MANEUVER  (Y.P.=1,W.L.T.=2,H.T.=3)  :" 

3930=  READ*,AA2 

3940=  PRINT*,"  " 

3950=  PRINT*,"  ENTER  INPUT  COMMAND" 

3960=  PRINT*,"  FINAL  VALUE,  RAMP  TIME  (UNIT  STEP  =  1,0)" 

3970=C 

3980=  DO  nil  I  =  1,L 

3990=  KVR(I)=0 

4000=  SL0PE(I)=0.0 

4010=  VST(I)=0.0 

4020=  KVST(I)=0 

4030=  PRINT*, I,"  :  " 

4040=  READ« ,VC( I ) ,VDT( I ) 

4050=  VC(I)=VC(I)/57.2957795 

4060=  IF<VDT(  I )  .E(3.0)  GO  TO  1160 

4070=  KVR(I>=1 

4080=  SLOPE(I)=VC(I)/VDT(I) 

4090=  GO  TO  nil 

4100=  1180  PRINT*,"  ENTER  STEP  RESPONSE  END  TIME  (INFINITY=0)  :" 

4110=  READ*,VST(I) 

4120=  IF(VST(I) .EQ.O)  GO  TO  1185 

4130=  KVST(I)=1 

4140=  1165  V(I)=VC<I) 

4150=  nil  CONTINUE 

«0=  TDT=VDT(1) 

W0=  IF<AA2.Ea.2)  TDT*VST<3) 

4180=  FV=VC(1 >*57.2957795 

4igo=c 

4200=  PRINT*,"  " 

4210=  PRINT*, "ENTER  FINAL  TIME  AND  TIME  STEP  I" 

4220=  READ*,TF,DELT 

4230=  PRINT*,"  " 

4240=  NEQN=NL 

4250=  T=0 

4280=  T0UT=0 

4270=C 

4280=  DO  1113  I=1,NL 

4290=  >{(I)=0 

4300=  Y(I)=0 

4310=  1113  U(I)=0 
4320=C 

4330=  DO  1112  1=1,5 

4340=  1112  ER(I)=0 
4350-C 

4380*  RR(1)=0 

4370=  RR(2)=0 

4380*  RELERR=lE-06 

4390=  ABSERR=lE-06 

4400=  IFLAG=1 

4410=  NSTEP*INT(TF/DELT)+1 
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4420=C 

4430= 

4440  = 

.^450= 

>460  = 

131 

4470= 

132 

4480= 

4490  = 

4500= 

125 

4510= 

4520  = 

4530  = 

126 

4540= 

4550  = 

127 

4560=C 

4570=C 

4580= 

4590=C 

4600= 

4610=C 

4620= 

4630= 

4640  = 

4650  = 

4660= 

4670= 

1155 

4680= 

4690= 

1150 

4700 =C 

(T710= 

4720=C 

4730= 

4740=C 

4750= 

4760= 

1120 

4770=C 

4780=C 

4790=C 

4800  = 

4810= 

1136 

4820=C 

4830=C 

4840= 

4850* 

4860  = 

4870=C 

4880= 

133 

4890= 

1128 

4900=C 

4910= 

4920  = 

4930= 

135 

4940=C 

4950= 

.660= 

136 

IF(AAl.EQ.l)  GO  TO  131 


PRINT*,"  T  OUTPUTS" 

GO  TO  132 

PRINT*,"  T  CONTROLS" 

PRINT*,  " - " 


IF(AAl.EQ.l)  GO  TO  125 
IF(AA2.EQ.l)  GO  TO  126 
PRINT  9980, TOUT, (U<I), 1=1, L) 

MRITE(4,9980>  TOUT, (U( 1 ), 1=1 ,L) 

GO  TO  127 

PRINT  9981,T0UT, (Y( I) ,I=1,L) ,RR{ 1) ,RR(2) 
WRITE (4, 9981 )  TOUT, ( Y ( I ) , I = 1 , L) , RR ( 1 ) , RR( 2 ) 
CONTINUE 


DO  1115  K=1,NSTEP 
TOUT=TOUT+DELT 
DO  1150  1  =  1, L 

IF(KVR( I ) .EQ.O)  GO  TO  1155 
0(I)=SL0PE(I)*T0UT 
IF<TOUT.GT.ODT<I))  g(I)=OC(I> 

GO  TO  1150 

IF(KOST< I ) .EQ.O)  GO  TO  1150 
IF<TOUT.GT.VST(I) )  0(I)=0.0 
CONTINUE 

CALL  ODE ( F , NEQN , X , T , TOUT , RELERR , ABSERR , I FLAG , WORK , I WORK ) 

IF(AA1.EQ.2)  GO  TO  133 

DO  1120  1  =  1, L 
U( I )=X( I3+I ) 


DO  1136  1  =  1, L 
UC(I)=U(I)*57. 29578 


PRINT  9980, TOUT, (UC<I), 1=1, L> 
WRITE(4,9980)  TOUT, (UC( I ) , 1=1 ,L) 
GO  TO  134 

DO  1128  1  =  1, L 
Y(I)»0 

DO  135  1  =  1, L 
DO  135  J=1,NL 
Y<I)=Y(I)+CEE(I, J)*X< J) 

DO  136  1=1, L 
Y<I)=Y(I)*57. 2957795 
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4970=( 

4980» 

4990= 

.  ■  ”000  = 
c>010= 
5020=( 
5030= 
5040  = 
5050= 
5060= 
5070= 
5080= 
5090= 
5100  = 
5110* 
5120»( 
5130  = 
5140  = 
5150* 
5160  = 
5170  = 
5180= 
5190= 
5200= 
5210= 
5220- 
5230= 
5240= 
^^50=C 
^260= 
5270= 
5280= 
5290= 
5300  = 
5310  = 
5320  = 
5330=( 
5340= 
5350= 
5360= 
5370  = 
5380= 
5390= 
5400= 
5410=( 
5420  = 
5430= 
5440- 
5450= 
5480= 
5470= 
5480= 
5490= 
5500- 
^TIO-C 
j520» 


134 

1115 

9980 

9981 


KK=2 

IF(TOUT.EQ.TDT)  KK=1 
IF<AA2.EQ.l)  GO  TO  137 
IF(AA2.EQ.2>  GO  TO  139 

IF(TOUT.GT.TDT)  GO  TO  141 
ER( 1 )=ER( 1 )+Y< 1 )*DELT*KK 
ER(2)=0 
ER(4>=0 

ER ( 3 ) =ER ( 3 ) +ABS ( Y ( 2 ) *DELT*KK ) 

E4= ( Y ( 1 ) -FV ) *DELT*KK 
ER(4)=ER(4)+ABS(E4) 

ER ( 5 ) =ER ( 5 ) +AB3 ( Y ( 3 ) #DELT*KK ) 

GO  TO  142 

IF(TOUT.GT.TDT)  GO  TO  138 
ER ( 1 ) =ER  < 1 ) +Y ( 1 ) *DELT*KK 
ER ( 2 ) =ER  <  2 ) +Y  <  3 ) *DELT»KK 
ER<4)=0 
ER(5)=0 

ER ( 3 ) =ER ( 3 ) +AB3 ( Y ( 2 ) *DELT*KK ) 

E4-<Y( 1 )-FO)*DELT*KK 
ER(4)=ER(4)-i-ABS<E4) 

ER  <  5 ) =ER ( 5 ) +ABS ( Y ( 3 ) *DELT*KK ) 

RR ( 1 ) »RR ( 1 ) +Y ( 3 ) *DELT 
RR(2)=RR( 1 )+Y( 1) 

GO  TO  142 

IF<TOUT.GT.TDT)  GO  TO  140 
ER< 1 )=0 

ER ( 2 ) *ER ( 2 ) +Y ( 3 ) *DELT*KK 
ER<5)=0 

ER ( 3 ) =ER  <  3 ) +ABS ( Y ( 2 ) *DELT»KK ) 
ER(4)=ER(4)+ABS<Y( 1 )#DELT«KK) 

ER ( 5 ) =ER ( 5 ) +ABS ( Y ( 3 ) *DELT*KK ) 

IF<AA2.Ea.l)  GO  TO  128 
PRINT  9980rT0UT,  (Yd)  rI»l,L) 

WRITE(4,9980)  TOUT, < Y( I ) r 1=1 »L) 

GO  TO  129 

PRINT  9981 , TOUT,  (Yd)  ,  1  =  1  ,L)  ,RR(  1 )  ,RR(2) 
WRITE (4, 9981)  TOUT,  (Yd  )  ,  1  =  1  ,L)  ,RR(  1  )  ,RR(2) 
CONTINUE 

CONTINUE 

CONTINUE 

F0RMAT(5X,F5.2,3(3X,F10.4) )  . 

FORMAT ( 5X , F5 . 2 , 3 ( 3X , F 1 0 . 4 ) , 2X , F7 . 4 , 2X , F7 . 4 ) 

PRINT*, " - 

PRINT*,"  " 

IF(AAl.EQ.l)  GO  TO  146 
DO  147  1-1,5 
ER(I)=ER(I)/KK 

PRINT*,"  BETA  COMMAND  =",ER(1) 
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ER(2) 


5530= 

PRINT*, " 

R  COMMAND 

5540= 

PRINT*, " 

PHI  =",ER( 

5550= 

PRINT*, " 

BETA  TRANS 

;*5S0= 

PRINT*," 

R  TRANS.  = 

3570= 

PRINT*, " 

II 

5580= 

146  CONTINUE 

5590= 

GOTO  4710 

5600=C 

seiosc 
5620sC 
5630= 
5640  = 
5650= 
5660= 
5670= 
5680= 
56S0=C 
5700=C 
5710=C 
5720= 
5730= 
5740= 
5750= 
5760= 
5770= 
5780= 
5790=C 
5800= 
^810= 
Ofa20= 

5830= 
5840  = 
5850  = 
5860  = 
5870  = 
5880=C 
5890  = 
5900  = 
5910= 
5920= 
5930= 
5940= 
5950  = 
5960= 
5970=C 
5980=C 
5990  = 
6000= 
6010= 
6020= 
6030= 
6040= 
6050= 
6060=C 
.  .  :  070«C 
■w080= 


9999 


CONTINUE 
PR I NT* r"  " 

PRINT*, "ALL  INFORMATION  IN  PORTER  HAS  BEEN  SAVED  IN  LOCAL" 
PRINT*, "FILE  -  PLOT" 

STOP 

END 


SUBROUTINE  MADD ( A , B , L ,N , C ) 

DIMENSION  A( 10, 10) ,B( 10,10) ,C( 10, 10) 
DO  100  1=1, L 
DO  100  J=1,N 
C<I, J)=A<I,J)+B(I, J) 

RETURN 

END 

SUBROUTINE  READM<L,N,A) 

REAL  A<10,10) 

PRINT*,"  " 

DO  100  1=1, L 
PRINT*,"  ROW  ",I," 
READ*,<A(I,J),J=1,N) 

RETURN 

END 

SUBROUTINE  PRINTM ( L ,N , A ) 

REAL  A(10,10) 

PRINT*,"  " 

DO  100  1=1, L 

PRINT  200, (A( I , J) , J=1 ,N) 

F0RMAT(1X,8<F12.4,3X) ) 

RETURN 

END 


SUBROUTINE  MSUB(A,B,C,L,N) 

REAL  A( 10, 10) ,B( 10, 10) ,C< 10, 10) 
DO  100  I=1,L 
DO  100  J=1,N 
C(I,J)=A(I,J)-B(I, J) 

RETURN 

END 


SUBROUTINE  DSGNl (N,L, A,B,C,D,E,F) 
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6090=  REAL  A( 10, 10) ,B( 10, 10) rC( 10, 10) ,D( lOr 10) ,E( 10, 10) rF( 10, 10) 

6100=  N=4 

6110=  L=3 

120=  A(1,1)=0. 

13130=  B<1,1)=0.0  *B<1,2)=1.0  *B ( 1 , 3 )=. 03184253 

6140=  C( 1,1) =.04885258  »C(2,1)=0.0  *C(3,1)=0. 

6150=  D( 1 , 1 )=-. 37442438  $D ( 1 , 2 )=. 03182685  *D ( 1 , 3 )=-. 99948763 

6160=  D ( 2, 1 ) =-39.92833732  *D(2, 2) =-2. 1 1380581  $D ( 2 , 3 )=. 07969577 

6170=  D(3, 1 )=6. 23794780  $0(3,2)=-. 00631591  *0(3,3) =-.70740917 

6180=  00  100  1=1,3 

6190=  100  E(1,I)=0.0 

6200=  F(  1,1)  =*.053091 16  »F  ( 1 ,2)  =  .00463283  *F  ( 1 , 3 )  =  . 02501018 

6210=  F ( 2, 1 ) =9. 79285579  ♦F ( 2 , 2 ) =-55 . 73875487  *F ( 2 , 3 ) =2 . 60326273 

6220=  F(3, 1 )=-5. 05439681  tF ( 3 , 2 ) =-l . 89267784  *F ( 3 , 3 ) =4 . 38800282 

6230=  PRINT*,"  " 

6240=  PRINT*,"  THE  All  MATRIX  16  -" 

6250=  CALL  PRINTM( 1 , 1 , A) 

6260=  PRINT*,"  " 

6270=  PRINT*,"  THE  A12  MATRIX  16  -" 

6280=  CALL  PRINTM ( 1 , 3 , B ) 

6290=  PRINT*,"  " 

6300=  PRINT*,"  THE  A21  MATRIX  IS  -" 

6310=  CALL  PRINTM(3, 1 ,C) 

6320=  PRINT*,"  " 

6330=  PRINT*,"  THE  A22  MATRIX  IS  -" 

6340=  CALL  PRINTM ( 3 , 3 , D ) 

6350=  PRINT*,"  " 

6360=  PRINT*,"  THE  B1  MATRIX  IS  A  1  BY  3  ZERO  MATRIX  " 

1^370=  PRINT*,"  " 

>^380=  PRINT*,"  THE  B2  MATRIX  IS  -" 

6390=  CALL  PRINTM ( 3 , 3 , F ) 

6400=  RETURN 

6410=  END 

6420=C 
S430=C 

6440=  SUBROUTINE  F(T,X,XP) 

6450=  REAL  X ( 10 ) , XP ( 10 ) , AY ( 10 , 10 ) , BEE ( 10 , 10 ) , 0 ( 10 ) 

6460=  COMMON  AY , BEE , 0 , NL , L 

B470=C 

6480=  DO  100  I=1,NL 

6490=  100  XP(I)=0 

S500=C 

6510=  DO  200  1=1, NL 

6520=  DO  200  J=1,NL 

6530=  200  XP(I)=XP(I)+AY(I,J)*X(J) 

6540=C 

6550=  DO  300  I=1,NL 

6560=  DO  300  J=1,L 

6570=  300  XP(I)=XP(I)+BEE(I,J)*M(J) 

6580=C 

6590=  RETURN 

6600=  END 

6610=C 

6620=C 

-■330»*E0R 

6640=*E0F 
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3300= 

3999  PRINT*,"  " 

3310= 

PRINT*, "PRINT 

C-L 

MATRICES  (1=YES)  I" 

3320  = 

READ*,X25 

-30  = 

PRINT*,"  " 

w-40= 

PRINT*, "ENTER 

THE 

GAIN  CONSTANT  ,  G  I" 

3350  = 

REAO*,G 

33B0=C 
3370=C 
3380=C 
3390=  9977 
3400= 

3410=  4000 
3420=C 
3430  = 

3440  = 

3450= 

3460  = 
3470=C 
3490= 
3490=C 
3500  = 

3510  = 

3520=  4090 
3530=C 
3540  = 

3550=  4095 
3560=C 
3570=  4100 
2Sfl0=C 


3600= 
3610= 
3620  = 
3630  = 
3640= 
3650= 
3660= 
3670= 
3680  = 
3690= 
3700  = 
3710= 
3720= 
3730  = 
3740  = 
3750= 
3760= 
3770= 
3780= 
3790  = 
3800= 
3810= 
3820= 
3830  = 
r‘  40= 


4190 


4195 


4200 


=  4290 


FORM  CL08ED  LOOP  A-MATRIX  (AY) 

DO  4000  1=1 fN 
DO  4000  J=1,L 
AY(I, J)=0. 

CALL  0MULFF(B2fKl >LrL,L,8r8fB2Kl rBr lER) 

CALL  VMULFF(B2»K0,L,L,Lr8r8,B2K0,8, lER) 

CALL  VMULFF(B2K0,F1 , L , L , NM r 8 , 8 , B2K0F1 , 8 , lER ) 
CALL  VMULFF<B2K0,F2,L,L,Lf8,8,B2K0F2,8, lER) 

00  4100  1=1, L 

DO  4090  J=1 ,L 

BKF2 (I , J ) =G*B2K0F2 ( I , J ) 

BKl(IrJ)=G«B2Kl(I,J) 

DO  4095  KA=1,NM 

BKFl ( I ,KA)=G«B2K0F1 ( I ,KA) 

CONTINUE 

CALL  MSUB(A21,BKF1,THRTW0,L,NM> 

CALL  M8UB(A22,BKF2,THRTHR,L,L) 

DO  4200  I=1,L 

DO  4190  J=1,L 

KA=I+N 

LA=J+N 

AY(KA, J)=BK1(I, J) 

AY(I,LA)=-F2(I, J) 

AY(KA,LA)=THRTHR( I , J) 

DO  4195  J=1,NM 
LB=J+L 

AY(I,LB)=-F1(I, J) 

AY<KA,LB)=THRTWO< I , J) 

CONTINUE 

DO  4300  I=1,NM 

DO  4290  J=1,L 

LA-I+L 

KA»J+N 

AY(LA,KA)=A12(I, J) 
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3850= 

3860  = 

3''70=  4285 
;:;:JJ0=C 
3880=  4300 
3800=C 
3910  = 

3920= 

3930  = 

3940  = 

3950= 

3960=  4310 

3970=C 

3980=C 

3990=C 

4000= 

4010  = 

4020=  4400 

4030=C 

4040= 

4050=  4410 

4060=C 

4070= 

4080= 

4080= 

4100=  4420 
4110=C 


4150= 


4160=  4510 
4170=C 
4180=C 
4190=C 
4200  = 

4210= 

4220=  4500 
4230=C 
4240= 
4250=C 
4260  = 

4270=  4520 
4280=C 
4290  = 

4300=  4530 
4310=C 
4320=  4515 
4330=C 
4340= 

4350= 

4360= 

4370  = 
-‘•^80=C 


DO  4295  J=1,NM 
KB=J+L 

AY(LArKB)=All(I, J) 
CONTINUE 


NL=N+L 

IF  (X25.NE.1)  GOTO  4310 
PRINT*,"  " 

PRINT*, "THE  CLOSED  LOOP  AY  MATRIX  IS  - 
CALL  PRINTM(NL,NL,AY) 

CONTINUE 

FORM  CLOSED  LOOP  B-MATRIX  (BEE) 

DO  4400  I=1,NL 
DO  4400  J=1,L 
BEEd,  J)=0.0 

DO  4410  1  =  1, L 
BEEd  ,  I )  =  1 .0 

DO  4420  I=1,L 
DO  4420  J=1,L 
GB2K0  d , J ) =G*B2K0  d , J ) 

BEE<N*1  ,  J)=GB2K0d  ,  J) 

IF  <X25.NE.l)  GOTO  4510 
PRINT*,"  " 

PRINT*, "THE  CLOSED  LOOP  BEE  MATRIX  IS 
CALL  PRINTM(NL,L,BEE) 

CONTINUE 

FORM  CLOSED  LOOP  C-MATRIX  (CEE) 

DO  4500  I=1,L 
DO  4500  J=1,NL 
CEE( I , J)=0.0 

DO  4515  1  =  1, L 

DO  4520  J=1,NM 
CEEd,L+J)=Cld,  J) 

DO  4530  J=1,L 
CEEd  ,N+J)=C2d  ,  J) 

CONTINUE 

PRINT*,"  " 

IF  (X25.NE.1)  GOTO  4610 

PRINT*, "THE  CLOSED  LOOP  CEE  MATRIX  IS 

CALL  PRINTM(L,NL,CEE) 
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